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1. INTRODUCTION AND SUMMARY 


The General Interpolants Method (GIM) code was developed to analyze 
complex flow fields which defy solution by simple methods. The code uses 
numerical difference techniques to solve the full three-dimensional time- 
averaged elliptic Navier -Stokes equations in arbitrary geometric domains. 

The equations are cast in strong Conservation law form and written in an 
orthogonal Cartesian coordinate system. Included are a continuity equation 
for global mass conservation, three components of momentum conservation, 
total energy conservation and an equation for conservation of individual specie 
of a binary gas. Pressure is related to the conservation variables through 
the ideal gas law for a binary mixture. A generalized geometry package is 
used to model the flow domain, generate the numerical grid of discrete points 
and to compute the local transformation metrics. Computation is done in 
physical space by explicit finite -difference operators. The GIM approach 
essentially combines the finite element geometric point of departure with 
finite difference explicit computation analogs. This provides a capability 
which takes advantage of the geometric flexibility of an element description 
and the superior computation speed of difference representations. 

The numerical analogs of the differential equations are derived by 
representing each flow variable with general interpolation functions. The 
point of departure then requires that a weighted integral of interpolants be 
zero over the flow domain. By choosing the weight functions to be the inter- 
polants themselves, the GIM formulation produces identically the classical 
implicit finite element discrete equations. These forms are not used in 
the GIM code due to their fully implicit nature and inherent inefficiencies. 
Rather, the weight functions are chosen to be orthogonal to the interpolant 
functions which produces explicit finite difference type discrete analogs. 

By appropriate choice of constants in the weight functions, the GIM be- 
comes analogous to such finite difference schemes as centered, backward, 


forward, windward and multi-step predictor-corrector schemes such as the 
MacCormack method. The GIM analogs, however, are automatically produced 
for arbitrary geometric flow domains and; hence is a more general point of 
departure and provides greater flexibility in choosing difference schemes. 

A motivation for developing this code on these; principles was to pro- - 
vide an analytical tool.which is more user oriented than the barsic i;esearch 
tools which exist., A fully production-line code to solve the complex Navier- 
Stokes equations does not exist today. In developing the GIM code, an attempt 
was made to bridge the gap somewhat between the p.ure research codes and 
the, ultimate production tool. The code was originally developed for a GDC 
7600 computer system. It has been subsequently converted to vector FORTRAN 
for the CDC STAR- 100 system at NASA Langley Research Center. Reference 1 
provides documentation for the GIM/STAR code designated version SE-1 (STAR- 
Elliptic No. 1). This version of the code has been used to compute a number 
of complex flow fields including nozzle flows for both subsonic and supersonic 
regimes, and two and three-dimensional Sc ramjet exhaust flow simulations 
(Ref. 2). 

The current contract work involves utilization and extensions of the 
GIM/ STAR code. The objectives of the study are to: 

• Compute flow fields in supersonic inlet configurations 
using the SE-1 code 

• Upgrade the technical capability of the SE-1 code 

• Develop a hyperbolic and a parabolic version of the GIM/ 

STAR code to supplement the elliptic capability. 

This report presents the progress to date on the development and application 
of the GIM/STAR code. 

Section 2 presents the results of an application of the code to a two- 
dimensional supersonic inlet. The calculation was started upstream of the 
compressioh surface which turns at 25 deg to the horizontal. The Mach = 5 
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freestream flow generates a bow shock off the leading edge of the ramp. The 
calculation involved two primary considerations; (1) determine the amount 
of flow captured and the amount spilled into the freestream and (2) compute 
the inlet flow field and predict the shock wave/boundary layer interaction. 

The problem was run in two parts with the GIM code on the STAR machine. 

The ingested flow was determined (inyiscidly) first and found to be 66% of 
the incoming stream. This agrees well with the numbers for which the simu- 
lated inlet was designed. The flowfield distribution at the nozzle entrance 
was then used to drive the internal flow allowing the performance parameters 
to be determined. The flow angularity produces a shock wave off the cowl 
lip which propagates into the nozzle. The ultimate interaction of this shock 
and the laminar boundary layer on the upper propulsion surface were com- 
puted. All shock waves were determined using the "capture" mode of calcu- 
lation. Section 2 shows the computed solution for the spillage part of the flow 
and for the internal nozzle portion. The separation of the boundary layer due 
to the adverse pressure gradient is clear from the velocity and pressure 
contour plots. Radial distributions of the steady state flow field are given 
and a "time" history of the shock/boundary layer interaction calculation is 
also shown. 

Section 3 of this report describes an investigation of linearized block 
implicit (LBI) finite difference schemes for the GIM code. The current 
explicit MacCormack schemes are relatively efficient for flows with in- 
viscid boundary conditions. In anticipation of other requirements to com- 
pute three-dimensional viscous flows, the necessity of eliminating the explicit 
stability limit becomes apparent. However, the extreme inefficiencies in- 
herent in "fully" implicit methods, due to the large band-width matrices, 
make them unrealistic for large three-dimensional viscous flow problems. 

The most promising concept is the linearized block implicit (LBI), or approxi- 
mate factorization, schemes. These methods retain the Conservation Law 
equation form while "splitting" the spatial dependence in the manner of the 
ADI schemes. The resulting matrix bandwidth is once again small (usually 3) 
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and is practical to use. The study of LBI schemes in this work was con- 
centrated on: 

• Stability requirements of the block tridiagonal scheme 
of Beam-Warming (Ref. 3) 

• Accuracy of the L.BI scheme itself and more precisely, 
the accuracy and speed of linear equation solvers for 
vector machines 

• Shock wave resolution of LBI schemes used in a capture 
mode and artificial damping requirements 

• Techniques to vectorize LBI schemes for use on the 
STAR- 100 machine, 

'^The study was carried out with a one- dimensional code that uses the Beam- 
Warming formulation. 

Results of the LBI investigation are discussed in detail in Section 3, 

The stability of the scheme was found to be strongly coupled to the accuracy 
of the linear equation solver used and to the artificial damping added to the 
explicit side of the scheme. The "unconditional” stability indicated by the 
theory could not be achieved numerically using centered differences. Schemes 
based on one-sided windward differences did prove to be unconditionally stable. 
The LBI scheme was shown to be as good as the explicit MacCormack for reso- 
lution of shock waves. The overall conclusion of this part of the study is that 
LBI schemes appear to be very promising for three-dimensional viscous flows 
but they are not as outstanding as the literature indicates. 

The third part of this study reported here is the development of hyper- 
bolic and parabolic methods to supplement the elliptic code. Section 4 describes 
the details of the work on the GIM maching algorithms and the current status 
of the code. The basic idea of the GIM code marching scheme is to combine 
the classical parabolized Navier-Stokes methods with a ”quasi-time” relaxa- 
tion. The term "quasi-parabolic” (QP) will be used to refer to this algorithm 
although the scheme applies equally well to hyperbolic, supersonic inviscid 
flows. The QP algorithm is contrasted to a fully elliptic method in that down- 
stream effects cannot be felt upstream and that a full flow domain need not be 
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stored for the QP scheme. The QP algorithm is also contrasted to classical 
parabolic methods in that mixed subsonic/supersonic flows do not produce a 
multiple "decode” root and that real-wall no- slip boundaries can be treated 
with the QP algorithm. The equations are the classical parabolized Navier- 
Stokes but with a psuedo-time derivative added back to them. The solution is 
known at upstream data planes 1, 2, . . . .N-1 and the solution is sought at plane 
N with no influence from plane N+1. Time relaxation is used to solve for plane 
N from only the (converged) solution at upstream planes. Backward differ- 
ences (second order) are used, of course, in the quasi-marching coordinate. 

As the algorithm is formulated, either explicit or linear block implicit time 
relaxation can be incorporated. 

The resulting algorithm then requires much less computer storage than 
a GIM elliptic flow field calculation and does not have the "singularities" in- 
herent in classical parabolic marching algorithms. The QP scheme has been 
coded and partially checked out on the STAR system. At the time of this 
writing, the GEOMETRY, MATRIX and INTEG modules of the SP-1 GIM code 
(STAR Parabolic, Version 1) have been run successfully for several sample' 
cases . 


Some details of the current contract work are appended. The most cur- 
rent version of the GIM elliptic code (SE-2) is discussed in Appendix A authored 
by L. W. Spradley. Differences in SE-1 and SE-2 are described and reasons for 
the changes explained. New INPUT data sheets for SE-2 are given to replace 
the ones in the "Blue Book" (Ref. 1). This basic guide should be used in con- 
junction with the Blue Book for inputting the GIM code on STAR. Appendix B 
by Jurgen Thoenes, contains a derivation and list of the complex linearization 
Jacobians for three-dimensional LBI schemes. The GIM-Marching code (SP-1) 
requires a special set of weight/shape functions. These are derived in Appendix 
C, which is authored by John F. Stalnaker. The final item to be covered here is 
a description of the vectorized linear algebraic equation solvers which were 
developed on the STAR system for use with the LBI schemes. The mathematical 
development and performance of several techniques, both direct and iterative 
are shown in Appendix D, authored by S, J. Robertson. 
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2. CALCULATION OF TWO-DIMENSIONAL INLET 
FLOWS WITH SPILLAGE 

2.1 INTRODUCTION 

Figure 2-1 shows the model two-dimensional supersonic inlet for 
which the flow field was computed using the elliptic GIM/STAR code. 

The compression surface makes a sharp 25 deg turn at x=0. It turns 
50 deg through a circular arc centered about x = 5 into the 25 deg ex- 
pansion surface. The expression surface and the lower cowl from the 
nozzle. The freestream flight conditions are also shown in Fig. 2-1. 

All flor variables are made dimensionless with the freestream quantities. 

For inlets with fixed geometry it is important to know the amount of 
flow captured by the inlet and the amount that spills into the freestream. Thus, 
special emphasis was placed on calculating the mass flow rate at the inlet 
throat (x = 5). The model inlet was designed inviscidly to capture 66.6% of 
the incident flow. 


It is felt that a brief history of the development and an outline of the 
pitfalls incurred in obtaining the final solution would be of benefit to future 
users of the GIM/STAR code. This discussion appears in Section 2.2. A 
complete analysis of the final solution is given in Section 2,3. These dis- 
cussions are divided into two parts: (1) the external flow field below the 

compression surface and including the freestream flow which spills below 
the cowl, and (2) the internal (nozzle) flow field. 


Use of trade names or names of manufacturers in 
this report does not constitute an official endorsement 
of such products of manufacturers, cither expressed or 
implied, by the National Aeronautics and Space Admin- 
istration. 
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Fig, 2-1 - The Model Two 




2.2 DEVELOPMENT OF THE SOLUTION 


2.2.1 External Flow Field 

To limit the problem size, the computational grid was constructed 
originally with the input boundary lying along the 35.7 degree bow shock line.. 

This resulted in computational difficulties with the grid points along a hori- 
zontal line from the leading edge of the cowl to the shock line. In order to 
wrap the grid around the cowl, a discontinuity in the shapes of the elements 
arose along this line. Relatively uniform rectangular elements were mated 
to severely skewed elements. It is believed that the computational problems 
arose from the finite difference analogs generated along this line of nodes. 

These improper influences were caused by either sharp discontinuities in 
the transformation metrics or inadvertent extrapolation in the transformations. 

As a result, the post- shock grid was abandoned and it was decided that 
the geometry should be constructed to allow the bow shock to be captured. The 
analogs for this grid were thoroughly examined using a coarsely spaced version 
of the final computational mesh (shown in Fig. 2-2). A "double-valued" node (L.e., 
two nodes at the same spatial location) was used to allow the proper splitting of 
the flow at the cowl Up, Due to the small shock angle at the bow which would not 
permit a sufficient number of nodes between the shock and the surface, the first 
eight nodes along the compression surface were held fixed at the inviscid post- 
shock conditions. This eliminated numerical disturbances which were otherwise 
generated at the bow and propagated downstream leading to instabilities. 

2.2.2 Internal Flow Field 

As originally modeled, the upper body of the inlet had a sharp 50 deg 
expansion at x= 5. In the initial inviscid analysis of the nozzle it was found 
that the flow overexpanded around this turn leading to pressure undershoot 
and instability. The sharp turn was rounded to alleviate this over expansion. 
However, subsequent analysis revealed the problem to be excessive damp- 
ing on the continuity equation. This resulted in an artificial dissipation 


9 



o 



Fig. 2-2 - Computational Grid for External Flow in the Two-Dimensional Inlet 


of mass away from the wall. Redaction of this damping allowed successful 
computation of the expansion; however, the rounded surface remained. The 
primary difficulty with the nozzle calculation was with the inviscid treatment 
of the expansion surface. This was first indicated by the failure of the in- 
viscid SEAGULL code (Ref, 4) to converge in the nozzle. Imposing a viscous 
boundary layer on the upper wall allowed the GIM code to develop a strong 
shock -boundary layer interaction which made evident the fallacy in the in- 
viscid treatment. 

In arriving at the final solution it has become increasingly clear that 
solutions with the GIM code are strongly dependent on two factors: (1) the 
structure of the computational mesh, and (2) proper modeling of the physics 
of the problem. 

2,3 RESULTS AND DISCUSSION 
2.3.1 External Flow Field 

The 3557 node computational grid for the external flow field is shown 
in Fig. 2-2. The solid boundaries were treated inviscidly. The USERIP option 
in the GIM/STAR code was used to initialize the flow field in order to lay in the 
bow shock as closely as possible to the inviscid 35.7 degree line. The solution 
converged to steady state in 900 iterations. The integrated mass flow rate 
indicated that the inlet captured 66.5% of the incident mass flow which com- 
pared almost exactly to the theoretical value. Figures 2-3 through 2-5 show 
the velocity vectors, pressure and Mach number contours for the complete 
flow field. Figure 2-6 shows a comparison of the mass flow rate (m = pu) 
across the inlet plane as calculated by the GIM/STAR code to that calculated 
by the inviscid SEAGULL code (Ref. 4) for a similar inlet with the same im- 
posed capture rate. The agreement is excellent with the only apparent differ- 
ences resulting from the different treatment of shocks in the two codes. For 
computational economy the deteched shock effects at the cowl lip were not 
treated here. Rather, after 100 iterations the values of the flow variables at 
the lip node were held fixed at attached post- shock conditions determined 
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Vertical Position 



Fig, 2-6 - Mass Flow per Unit Area vs Vertical Position Across the Inlet 
Plane — (Station 5) Comparison with SEAGULL Solution (Ref. 4), 


from conditions immediately upstream of the cowl lip. Further, when the 
lip shock was allowed to detach, the blunt body effects did not extend beyond 
one grid width from the Up. 

2«3.2 Internal Flow Field 

As noted in Section 2,2.2, it was necessary to impose a viscous boundary 

layer on the upper wall of the nozzle to obtain the solution. The boundary 

layer profile at the throat was estimated by a quadratic laminar profile with 

4 

a Reynolds number of 1 x 10 . The inlet flow variables were input by linear 
interpolation of the external flow results and the remaining nodes exterior 
to the boundary layer were initialized via the USERIP option by an isentropic- 
like area expansion along streamwise rows of nodes. It was found that allow- 
ing the code to develop the shock in the nozzle was preferable to estimating 
the shock position as was done in the external flow calculation. Figure 2-7 
shows the 3000 node computational mesh for the nozzle. The solution con- 
verged to steady state in 1200 iterations. The final velocity vectors, pres- 
sure and Mach contours are shown in Figs. 2-3 through 2-5. Figures 2-8 
through 2-20 show the time development of the solution from iteration 0 
through 1200. Figures 2-21 through 2-24 show variations of the Mach number 
and pressure in the nozzle compared with the available SEAGULL results. It 
is apparent from these last figures that the boundary layer is artifically too 
thick (a result of the choice of Reynolds number). However, the result ob- 
tained provides considerable insight into the physics of the problem as well 
as the reasons behind the failure of the Lnviscid analysis. 
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Fig, 2-7 - Computational Grid for Internal Flow in the 
Two-Dimensional Inlet. 
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on 600 ). 



Fig. 2-15 - Two-Dimensional Spillage Problem (Viscous Nozzle 

V = 4.87; Iteration 700). 
max 





Fig. 2-16 - Two-Dimensional Spillage Problem (Viscous Nozzle; 
V = 4,90; Iteration 800). 

XXXn.X 
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Fig. 2-17 - Two-Dimensional Spillage Problem (Viscous Nozzle 

V = 4. 91; Iteration 900). 
max 





Fig. 2-18 - Two-Dimensional Spillage Problem (Viscous Nozzle; 

V = 4. 91; Iteration 1000). 

max 
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Fig. 2-19 - Two-Dimensional Spillage Problem (Viscous Nozzle; 

V = 4. 91 ; Iteration 1100 )* 
max 
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Fig. 2-20 - Two-Dimensional Spillage Problem (Viscous Nozzle; 
V = 4 . 91 ; Iteration 1200 )* 



VertLcal Distance from Lower Wall (Cowl) 



Fig. 2-21 - Mach Namber Variations Across the Nozzle. 
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Fig, 2-23 - Pres 



e Variations Across the Nozzle. 
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Fig. 2-24 - Pressure Variation on Upper 
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3, INVESTIGATION OF LINEARIZED BLOCK 
IMPLICIT METHODS FOR THE GIM CODE 


3.1 INTRODUCTION 

Numerical solution of the unsteady Navier -Stokes equations by explicit 
finite difference techniques has a number of disadvantages. The most serious 
one, from a practical engineering viewpoint, is the small time steps which are 
usually required to maintain stability. Computation of boundary layer flows at 
high Reynolds number requires fine grids near solid boundaries, hence very small 
time steps and long computer run times. One apparent cure for these difficulties 
is the use of implic it methods some of which are unconditionally stable for any 
size time step. These schemes are not without problems of their own in terms 
of their practical use. Among the major difficulties are the following; 

1. Implicit finite differences, in general, lead to systems of nonlinear 
algebraic equations when applied to the Navier-Stokes equations. 

These must either be solved directly or linearized in some manner. 

2. Direct linearization, via classical ADI processes, will destroy the 
Conservation Law Form of the Navier-Stokes equations and hence 
shock capture algorithms cannot be used. 

3. Multi-dimensional implicit methods lead to very large systems of 
simultaneous algebraic equations. Even for linear systems, the 
efficient solution is not practical due to large size of the matrix 
coefficients. 

4. Fully implicit methods cannot be programmed for efficient use on 
advanced vectorized machines such as the STAR, ILLIAC, or NASF. 

Numerical treatment of the steady state parabolic form of the Navier- 
Stokes equations face many of the same difficulties as the elliptic form. The 
spatial marching step size is constrained by the small grid required to resolve 
boundary layers normal to a solid wall. Marching downstream great distances 
can result in impractically long run times. Implicit finite differences have 
the potential to eliminate the difficulties mentioned above. 
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3.2 ONE -DIMENSIONAL UNSTEADY DEVELOPMENT 


The first item to be developed Is the formulation of an implicit scheme 
which results in a linear algebraic system yet retains the conservation law 
form of the Navier -Stokes equations. This idea can be explored by consider- 
ing the equations in one space variable, x, and the time coordinate, t. 


At 


X 


Ax 


Direct linearization is usually done by "lagging" certain of the nonlinear con 
tributions by one time step. This destroys the conservation nature of the 
Navier-Stokes equations. 


The case considered here is an elliptic boundary value problem in space 


Governing Equations 


U = 


in time 

. The 

equations 

^ + 

0E 

9t 

3t + 

ax " 

0X 



/pv 

(pv) 

E = 

1 pv + P 

\P<?/ 


\(p<?+ P) 


( 1 ) 
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where 


p = mass density 
P = pressure 
jl = viscosity parameter 
(2 p" + X) 

t = time coordinate 


V = flow velocity 

^ = total energy 

k = thermal conductivity 

X = space coordinate 


P 


(Y-l)p[^-- vV2 


ideal gas law 


General Finite Difference Form 


This analysis will use the ’’delta" form of the flow variables 


AU^ = 


- 

AE*^ = 


- E^ 

At^ = 

n-H 

T 

n 

- T 


where n is the time step index. All data are assumed known at n= 0, Solving 
for AU^ then allows the data at level n to be advanced to level n+1; 


U 


n+1 


= + AU^ 


The class of finite-difference schemes considered can be written as follows; 


Au" = ^ (AU") + — (U 

i+€ at ^ ^ i+e at 


Au”' ^ + o|(0 - ” - e ) At^ + Af 


( 3 ) 


The parameters 0, c are used to generate a specific type of scheme. For 
0=0, the scheme is fully explicit; 9 > 0 gives an implicit method. If € = 0, 
the scheme requires two data sets of storage at time levels n, n+ 1. If 
€ > 0, then three levels are required to be stored, n - 1, n, n+1. 
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If 0 = 1/2 + e, the scheme Is second order accurate; and is first order 
otherwise. In this work, we are primarily concerned with 

0=1 e = 1/2 

which is a second order, implicit, three level scheme. 


Development of the Scheme 

The differential equation (1) is substituted into the general scheme (3) 
to get the following form: 


AU 


n _ 6 At 
1 + e 


a 

. A n A n. 

At 

9 , n n 

9x 

(-AE + At ) 

^ 1 + e 

^{-E +T)] 


+ 1 (4) 


By approximating the spatial derivatives, 9/0x, by finite-differences, we get 
a set of nonlinear algebraic equations. The incremental flow variables, AE^, 
At^^ are nonlinear functions (1) of the independent vector AU^. 

For our implicit scheme, this would require a simultaneous nonlinear 
algebraic equation solver. The best known methods are iterative ones which 
require long computer runs. 

For this work, we will perform a linearization as follows to obtain a 
Bot of linear algebraic equations and use matrix methods for their solution. 
The main idea here is to linearize the algebraic equations, but retain the 
fully conservative nature. 

Expanding E, t i-n a Taylor series, we get 

= e“ + (!§)" - u") + O(At^) 
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Ae” = j Au“ + O(At^) 

= a“ Au" + p(At^) 

and ^ 

= t“ + (1^) Au” + (^) AU" + 0(At") (5) 

or 

At* = P AU + R AU^ + 0(At ) 

where 

= 8U/0X 

The expression for At can be rewritten in a more convenient form by expand- 
ing the x-derivative to get 

At^ = (P - R^) AU^ O(At^) 

where 

R = 8R/8x 

X 


This form produces a linear system of equations with the same formal accu- 
racy (At^) as the nonlinear set. It does however, require evaluation of the 
Jacobians 


and 


A 


8E -p _ ^ 

au 9U 8U 

X 


8R/0X 


(6) 


39 



Putting the Taylor series (5) into the scheme (4) gives the following 
expression: 


AU 


n 


eAt a 

1+6 8x 

At 8 
1 + e 9x 


-a" Au" + (P - R^) Au" + (RAU)” 


(7) 


.n 


n. 


{-E + T ) + 


l + € 


AU 


n-1 


The last two terms on the right hand side of Eq. (7) are all explicit at time 
levels n, n-1. Denote this by D^, and write Eq. (7) as follows; 


au" + 


1 + e dx 


(A - P + R ) 
' x' 


n 


AU 


n 




1 = 


D 


n 


( 8 ) 


For convenience, let 


h = B = A - P + R 

1 + e X 


and write Eq. (8) as follows 


AU^ + h 


^ (b"au") - (r"au") 

9x 


= D 


n 


(9) 


To see that the form Eq. (9) may be useful, we will now write it for node point 
"i" in space and use second order centered finite differences 


dx 


I 


f. 


i+l 


- £. 


1-1 


2Ax 


+ 0( 


x^) 



^i+1 ■ ^ ^i ^i-1 
Ax" 


+ 0(Ax ) 


( 10 ) 
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With these difference expressions, Eq. (9) can be written as follows; 

r 


AU^^+ h 


^1+1 ^^i-n - ^1-1 

2 Ax 


Ax" 


= D 


n 


( 11 ) 


Combining coefficients of each AU^i AU^ji AU^_ ^ terms gives 


h ^n h _n \ . /t . 2h „n\ 

2Ax ^i+1 " A 2 ^i+l)^^i+l . 2 ^i ) ^^i 


Ax 


Ax 


( 12 ) 


h ^n h ^n \ ^n 

-T7T B- , 5 R. , ) AU. , = D, 

2Ax L-1 ^2 *•■'1/ ^“1 *■ 


(where I is the 3x3 identity matrix) 


Boundary values i= 1, and i = K must be treated separately due to the 
centered differences. For now we will let i = 2, 3, ... k - 1 and worry about 
boundary conditions later. 


The coefficients of the AU terms are 3x3 matrices which couple the 
three governing equations at each node point. There is an equation (12) for 
each node i = 2, 3, ... k - 1. 


To readily see the character of this system of linear algebraic equa- 
tions, let 


n 

-h 
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i+1 


(13) 
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pr in matrix notation: 



The system (14) will be termed " block tr idlagonal ." The individual 
matrices are full 3x3 arrays but they are arranged in a tridiagonal manner 
in the full matrix. The block arrangement occurs due to tlie linearization 
scheme used. This effectively couples the three differential equations at 
each node point. The boundary values for i= 1, and i = K have not been 
treated. This is an additional development item. 

The advantages of a system like Eq. (14) are: 

1. The Conservation Law form has been retained. 

2. Block tridiagonal systems are not much more costly to solve than 
pure tridiagonal systems. 

3. Operations like (14) can be vectorized for use on STAR -like ma- 
chines. 

Formulation of this scheme requires the analytical evaluation of the Jacobian 
matrices A, P, R. A brief look at these operations now follows. 


42 



Calculation of the Matrices 


The final matrices needed are L, M, N in Eq. (13). These are made up 
of combinations of B, R matrices from Eq. (5). 



,n ^ n 

» = A - P + R 

X 


A 


n 




R 


n 

X 



(15) 


The matrices A, R will have relatively simple elements (as we shall see), 
but the P, R^ matrices will be quite complex. For now we will assume that 
the viscous coefficients are constants; hence we will see that 


P - R =0 

X 

(See Beam-Warming paper. Ref. 3), 

We then need to analytically evaluate A, R, where 

BE. dT: 

au. ^ij du^ 

j 

The algebra for these operations is straightforward and is not included here. 
The final results are; 
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(16) 


where is the (constant) specific heat at constant volume, y is the ratio 
of specific heats and k is the thermal conductivity. 


Summary of Computational Procedure 


1. Set initial data at t = 0 for all nodes i= 1, 2, ... K. 

2. Form the vectors U, E, j. 

3. Compute D by explicit differences (Eq. (8)). 

4. Evaluate A, R matrices from Eq. (l6). 

5. Form the L, M, N matrices from A and R (Eq. (13)). 

6. Modify for boundary values. 

7. Call TRIDAG in the GIM code logic to solve the block 
tirdiagonal system for AU^. 

8. Advance solution vector to (n+1) 


U 


n+1 


= + AU^ 


9. Repeat the process to step 2 for a specified number of 
steps or until |AU^^1 < 6 for convergence to steady state. 
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3.3 ONE -DIMENSIONAL LINEARIZED BLOCK IMPLICIT RESULTS 


The procedure outlined in Sections 3.1 and 3.2 was subsequently coded 
and checked out. The equations were modified slightly to handle the problem 
of an expanding duct, quasi one-dimensional, by the inclusion of the area terms. 
This permits the computation of flows other than just the trivial case of con- 
stant property flow through a constant area duct. Three cases were con- 
sidered in order to check out and prove the method. Consider Fig. 3-1 where 
the simplest case is when the inflow conditions are fixed at the upstream end 
of the duct. For completely subsonic flow, elementary considerations indicate 
that the outflow at the downstream end of the duct has a unique solution. Con- 
sider for the moment that the flow is controlled entirely by the inflow conditions 
and the out-flow conditions are permitted to develop freely. Of course, it is 
known that physically one could change the back pressure at the downstream 
end and this would affect conditions at the upstream end. But, computationally, 
we specify the inflow conditions and therefore all the flow properties are 
uniquely determined. The same reasoning applies to the case where the flow 
is completely supersonic. In this case there exists the choking effect which 
means that when the back pressure is lowered below the limiting value no 
upstream effect is felt. If however, the back pressure is raised, the situa- 
tion develops where a normal shock moves into the duct with its positioning 
depending upon the back pressure. Thus for fixed inflow conditions for the 
supersonic case an unique solution depends upon the outflow pressure. Since 
so much is known analytically about this quasi one -dimensional case it was 
deemed a reasonable test with which to evaluate the linearized block implicit 
(LBI) scheme. 

In order that the LBI scheme could be applied to all three cases, i.e., 
including the strong shock case, pseudo viscous effects were included in the 
original coding in terms of numerical diffusion cancellation (NDC) terms. The 
first case computed was for fully supersonic flow through the expanding duct. 
The LBI scheme worked well reproducing the analytical results within about 
2% over the length of the duct. The case was initially run at a Courant num- 
ber of one. Subsequent runs were made at larger Courant numbers up to 3 
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Fig. 3- 1 - One Dimensional Duct Configuration 
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with good results. Increasing the CFL multiplier further caused rapid de- 
terioration of the solution and ultimate destruction of the case (it blows up). 

Theoretically the implicit solution should work for very large Courant 
numbers. Mathematically this is, of course, true but it ignores the physics 
of the situation. To verify that the solution was correct the complete deriva- 
tion was double checked, the coding was rechecked and nothing was found 
wrong. At first it was thought that the non-dominance of the main diagonal 
might be causing matrix ill- conditioning. The super- and sub-diagonals are 
both proportional to the step size while the main diagonal remains constant 
(at least for the inviscid equations). A natural conclusion might then be 
drawn that, as the step size is increased, non-dominance could occur such 
that the solution of the block matrices loses accuracy thus destroying the 
solution. 

To test out this theory some numerical experiments were carried out. 

First, an unblocked scalar matrix with three diagonals was used. A known 

solution was fed into the matrix reduction scheme and the non-dominance 

factors between the main and other diagonals were increased gradually. The 

case was run on the PDP-11, single precision arithmetic and inaccuracies did 

show up in the sixth place for even a Z to 1 non-dominance ratio. At 10^ to 1, 

9 

inaccuracies occur in the first and second places and at 10 to 1, order of mag- 
nitude inaccuracies were produced. Using double precision arithmetic on the 
PDP-11 or running the case on a GDC 7600 produced no inaccuracies whatsoever. 
Thus it is concluded that scalar matrices manipulated on high precision computing 
equipment have no accuracy problems associated with diagonal non-diminance. 

The same type of numerical experiments were then conducted with the 
block matrices. The CDC 7600 was used in order to eliminate any inaccu- 
racy due to less precise computing equipment. Non-dominance ratios ori the 
order of 10^ to 1 were necessary to generate errors in the fifth and sixth 
place. Since the suspected non-dominance caused by increasing the step 
size would only be of order 10, it is concluded that the reason the case would 
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not run at large Courant numbers Ls due to the problem physics and is not 
related to the mathematics of diagonal non-dominance. 

Subsequent consultations with NASA -Ames personnel (Robert Warming 
and Richard Beam) indicated that they saw no accuracy problems related to 
non -dominant diagonals and they believe the problem with using large Courant 
numbers is due to physically unrealistic propagation of pressure signals which 
then cause oscillatory behavior and eventually a negative pressure. Two 
different solvers were used to eliminate the possibility of an error in the 
coding. The two solvers, one from Lockheed-Huntsville and one from Ames 
Research Center, produced identical results, A fourth order damping term 
was appended to the RHS to help alleviate some of the oscillatory behavior. 

Ames indicated that in all their calculations with centered differences, 
fourth order damping was used. A fourth derivative term was therefore 
approximated and added to the RHS of the equations. The numerical diffusion 
cancellation terms were then dropped, except for the cases with shock. Use 
of the damping term eliminated some of the spatial oscillation but is highly 
dependent upon the value of an arbitrary coefficient which can vary between 
0 and 2. If too small a value is used the parameters oscillate, if too large 
a value is used the solution is overdamped and becomes linearized. A com- 
promise value used throughout this study was 0.1 which worked quite well for 
most of the cases analyzed. 

Another idea that was investigated involved the use of a MacCormack 
operator to compute the RHS. It is well known that the two-step MacCormack 
operator gives second order accuracy and is very stable. This scheme worked 
quite well and eliminated the necessity of including the fourth order damping 
term. 


To this point all the calculations were done using three point centered 
difference approximations to the derivatives. This results in the basic block 
tridiagonal scheme. One can also formulate the equation set based upon a 
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backward difference approximation which then results in a block bidiagonal 
scheme. This approach worked very well and, as is well known, has excellent 
stability characteristics. Its major limitation is that it is only first order 
accurate and generally is applicable only to supersonic flows. This scheme 
is inherently stable for any step size, and several cases were run at Courant 
numbers of 1000, 

Instead of using a pure centered scheme which has stability problems 
or a backward difference approximation that is only first order accurate, a 
combination of weighted differences was evaluated. Several combinations of 
weight factors were investigated, such as 2/3 centered plus 1/3 backward, 
and generally It was found that a slight increase in the Courant number could 
be obtained over that required for the pure centered scheme. Accuracy re- 
mained about the same as the centered differencing scheme. 

As Fig. 3-2 shows, the solution technique previously discussed produces 
very reasonable results including the location of the normal shock in the 
diverging duct. 
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FLg. 3-2 - Quasi-One Dimensional Results. 
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4. DEVELOPMENT OF GIM/STAR SPATIAL 
MARCHING ALGORITHMS 

4.1 INTRODUCTION 

The GIM/STAR SE-1 code treats the full elliptic flow field using explicit 
finite difference methods. This technique is applicable to a large range of 
fluid dynamics problems and has been successful in computing a number of 
these. The current code can be an "overkill” for some problems of interest 
in that a full elliptic treatment is not necessarily required. A parabolized, 
spatial marching algorithm could provide accurate flow fields for these situa- 
tions and would be considerably more economical. 


The elliptic code is constrained by two items which restrict its use on 
large three-dimensional viscous flows: 


1. The time step in explicit schemes is restricted by the 
CFL and viscous stability limit. This is usually controlled 
by the small grid sizes normal to no -slip boundaries. 

If inviscid, free slip boundaries can be used, i.e., ignore 
the boundary layer, then the severity of this constraint 
decreases. The implicit, linearized block methods 
described in Section 3 provide a possible remedy for 
the time step difficulty in the elliptic code, 

2. The large amount of data storage needed for three- 
dimensional viscous flows causes large "page faulting" 
on the STAR machine. Any finite difference method, 
explicit or implicit, still requires the large data base. 

A GIM/STAR code with a parabolic spatial marching 
algorithm would not attack as many kinds of problems 
as the elliptic version but would allow large three- 
dimensional viscous flows to be treated with no page 
faults on STAR. 


The intent of this research is to provide both an elliptic, time -dependent 
GIM/STAR code and a hyperbolic/parabolic spatial marching version. It 
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is not too difficult to conceive of the future codes which could contain switch- 
ing logic to automatically change from elliptic to parabolic etc., depending 
on the physics of the flow, however this will not be attempted here. The 
section will review the available classical parabolized methods and their 
problems, and then present an idea believed to be new for computing quasi- 
parabolic flows. 

4.2 CLASSICAL PARABOLIC APPROACHES 

Most of the literature on spatially marching schemes, hyperbolic or 
parabolic, treat equations which have been transformed to a Cartesian com- 
putation grid which is uniform. The space marching can then be done in 
much the same manner as time marching. This is a good approach if a 
single transformation exists for the full flow domain. The GIM code strategy 
has been to compute in the physical domain whereby completely arbitrary geom- 
etries can be treated. This approach presents a problem in developing a space 
marching algorithm, i.e., the fact that the geometry changes in the marching 
coordinate direction. This is akin to a GIM unsteady time marching scheme 
whereby the geometry is allowed to change with time. If we are to keep the 
GIM strategy of arbitrary geometries, then a space marching scheme must 
be developed which will account for the geometric variations in the stream- 
wise direction, i.e., non-uniform computational domain. 

The recent work of Roberts and Forester (Ref. 5) use a boundary -fitted 
computational mesh in a parabolic code for ducts of arbitrary cross section. 
Their algorithm for solving the equations appears to be a refinement of the 
classical method of Patankar and Spalding (Ref. 6). Rubin and Lin (Ref. 7) 
presented a nonlinear, iterative finite difference method for three-dimensional 
viscous flows. A parabolic method using a block implicit type scheme was 
given by Hirsh (Ref. 8). The solutions were restricted to supersonic flow 

(shear layers) of the free mixing type, Lubard and Helliwell (Ref, 9) calcu- 

\ 

lated flows on cone at angle of attack using a parabolized method. This paper 
discussed some of the inherent difficulties with singularities, ambiguities and 
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departure solutions which arise in parabolized algorithm. The paper dis- 
cusses explicit and implicit schemes for parabolic marching flows. 

Lin and Rubin (Ref. 10) presented a method using psuedo-time relaxation 
with a space-centered implicit differencing technique. They discuss many of 
the problems inherent in ’'pure” parabolic marching and show how time relaxa- 
tion can eliminate departure solutions. The GIM technique, although developed 
independently of Lin and Rubin, also employs time relaxation but with an ex- 
plicit, one-sided, predictor-corrector scheme and arbitrary three-dimensional 
geometries. The second order backward-forward, backward -backward explicit 
scheme of the GIM code is also a unique approach to parabolic marching solutions 

For problems in which viscous terms can be neglected entirely and the 
main flow direction remains supersonic, we would like the capability in the GIM 
code to resort to a simple hyperbolic algorithm. The classical methods pre- 
sented in the literature for parabolic and hyperbolic flows are drastically differ- 
ent because of the treatment of the pressure terms in the marching direction. 

As long as the flow is inviscid and supersonic, the axial pressure terms can be 
treated exactly. However, for subsonic flow, for example, several problems 
arise in applying a hyperbolic algorithm to the parabolic equations. This is, 
however, the approach that would be most general, if the "parabolic pressure" 
problem can be treated. 

Certain assumptions must be made in using a spatial marching technique. 

• There must exist a dominant flow direction in which to march. 

There can be no flow back upstream, i.e., no recirculation in the 
streamwise direction. 

• Stress terms are not allowed to act on the cross planes: i.e., there 
can be no second order terms (diffusion, viscosity) in the marching 
coordinate. 

• The downstream pressure field must not be allowed to propagate 
upstream. 

There are a number of strong implications in these assumptions. A super- 
sonic, inviscid flow satisfies them all. A supersonic viscous flow will con- 
form to the assumptions if the viscous terms are dropped in the marching 
coordinate direction. 
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Consider now the problem of spatial marching in a subsonic viscous 
flow. The first two of the assumptions can be met by simply not allowing any 
flow reversal problems to be attempted and dropping all streamwise diffusion 
terms. The downstream pressure field can still feed back through a subsonic 
stream. One obvious approach is to drop the streamwise pressure gradient 
term. This would satisfy the third assumption, but it appears a serious 
matter to simply drop this important term. 

Another approach commonly used is to provide a separate, explicit equa- 
tion for the pressure and use windward, one-sided differences. The most exact 
way is to compute the conserved flux parameters and then "decode" for the 
velocity, density and energy and compute the pressure from a state relation. 

The ideal gas law, a set of equilibrium thermodynamic relations or Boussinesq 
equations, is used to couple the state variables. Each of these approaches con- 
tains inherent difficulties which render their general use questionable. The 
following is a summary of some of these problem areas with classical parabolic 
schemes. 

Zero Axial Pressure Gradient 


This does not cause any significant numerical problems in computing a 
flow field. It does however create a major problem in that the computed answers 
are probably wrong for most flow fields. A mixed supersonic/subsonic flow, for 
example, with a shock wave crossing the flow field cannot be computed at all be- 
cause of the large axial (and radial) gradients. Some researchers still proceed 
to use this approach and try to justify it. 

Exact Pressure Treatment 


The rigorous way to compute the parabolized equations is to include the 
pressure in the conservation variable state vector for the momentum equations. 

A state equation can then be used to "decode" for the pressure. The advantages 
of the approach are that: (1) fully conservative differencing can be used; (2) shock 
capture algorithms are applicable; and (3) an auxiliary differential equation for 
the pressure is not needed. However, there are major problems with the "exact" 
treatment of the parabolic pressure. 
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• One-sided upstream differences must be used 

• The "decode" is ambiguous at Mach := 1 since two 
roots appear for the velocity (or pressure) 

• Real viscous, no- slip walls cannot be treated since 
the decode is singular, 

• Flows with a quiescent part, such as jets exhausting 
into an ambient, motionless atmosphere cannot be 
treated because of the singularity for zero velocity. 

Consider the two-dimensional parabolic system 


|E ^ ^ = 0 

a X oy 


where x is the marching coordinate, y the cross plane (or radial) coordinate, 
E is the state vector of conservation variables and F is a nonlinear (viscous 
plus convective terms) function of E. 


A typical state vector E for the parabolized Navier -Stokes equations is 


pu 

2 

pu + P 
puv 

{pS + P) u 


P 


O' - 1)P 


- 


2 2 
u +v 


Here, u is the axial velocity, v is the cross plane coordinate velocity, p is 
density, ^ is the total energy per unit volume and P is the pressure. Suppose 
a calculated value for the E vector exist at a plane X = X^, It is now required 
to "decode" for the primitive variables. The following is one decode pro- 
cedure that can be used in computer codes. 
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( 1 ) 


V = 


( 2 ) 

(3) 

(4) 

(5) 



p = Ej/u 


E 




p= (y=i)P 


e- 


2 ^ 2 
U + V 


Two problems are immediately obvious; 


• The radical in the u velocity decode causes an ambiguity. 
It can be easily shown that the correct decode is to take 
the + sign for u supersonic and the - sign if u is sub- 
sonic. In mixed flows, the sonic nature of a grid point 

is not known a priori. This Mach= 1 ambiguity prohibits 
a general parabolic marcher from being developed using 
the classical notions. See Section 5 of Appendix B. 

• The axial velocity, u, cannot be zero or the decode is 
singular. The axial component must be zero, however, 
if a real wall is to be put into the problem. All classical 
parabolic codes simply use some wall functions or resort 
to inviscid slip conditions to avoid the singularity. 


A third difficulty, which is not so obvious, is that attempts to use implicit 
methods to march the solution downstream often fail. The reason is that 
the boundary conditions are not treated exactly, and these errors build up as 
the streamwise coordinate is traversed. Often, the explicit differencing of 
points near the wall is used as a patchwork way of circumventing the boundary 
condition difficulty. 
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Explicit Treatment of Pressure 


This is the most widely used of the parabolic procedures and its origi- 
nation is usually attributed to D. B, Spalding. The idea is to provide an explicit 
differential equation for the pressure field in addition to the basic conservation 
laws. This is usually a Poisson-type relation obtained from combining con- 
tinuity and momentum equations. Satisfaction of local mass conservation is 
generally the criteria used for convergence of the elliptic Poisson equation. 

In general, a state vector will have the following appearance; 


E 


pu 

pu^ + P 
pu V 

(p<?+P)u 

P 


where the E^ component now represents the differential equation solution for 
pressure from whatever means. 

Now note the difference in the "decode" from the exact treatment case: 



V = 


p 



s 



- P/p 
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The Mach = 1 ambiguity is no longer present as the radical does not appear. 
Thus mixed subsonic/supersonic flows. can be coimputed without a priori knowl- 
edge of the Mach number. Note hpw.e.ver, that the. decode still contains the 
axial velocity in the denominator. Real solid walls cannot enter if viscous 
boundary conditions are used. 

This ” explicit pressure" treatment requires solution, at each plane, of 
a Poisson-type equation. Thus, an iteration between planes is required before 
movir^ on down to the next plane. Even with its inherent bad points, this ap- 
proach remains the most successful and widely used parabolic algorithm, 

4.3 THE QUASI-PARABOLIC IDEA 

The results of the initial investigation of a parabolic/hyperbolic GIM 
code led to the conclusion that there just is not a good approach being used 
today that fits the GIM code strategy. Three basic requirements were placed 
on a GIM/ parabolic algorithm: 

• The geometric treatment must be applicable to arbitrary 
shapes. 

• The same basic algorithm should be applied to both hyper- 
bolic and parabolic flows and be capable of eventual coupling 
with an automated algorithm for switching back and forth to 
the elliptic solver. 

• The algorithm should be readily vectorizable to realize the 
speed gain from using the STAR computer. 

In terms of a "classical" parabolized spatial marching algorithm, several 
geometric approaches were investigated. 

The first approach considered would generate the geometry plane by 
plane as the solution evolves, assembling the elements locally at each step. 

This would of course mean that the GEOMETRY module would be called at 
each integration step, thus coupling the geometry and the flow. An advantage 
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of this approach is that only the amount of geometry needed would be com- 
puted and stored at any cycle of the calculation. This would reduce the 
computer storage and reduce the large input/output problems. However, 
this approach would also require considerable reprogramming to make the 
GEOMETRY module of GIM a subprogram to the INTEGRATION module. 

A second approach appears to be a treatment of the geometry uncoupled 
from the flow. This means that all geometry, transformations and element 
assembly would be done before the flow field is integrated. The matrix data 
would be read from a stored file for each cross plane as needed. The ad- 
vantages of this approach are the geometry is computed only once for a 
given configuration, the geometry module can be separate (as it is now), 
from the integration module and the grid could be inspected for desirable 
character prior to computing an expensive flow field. Disadvantages are 
that the basic character of the flow must be analyzed a priori to place grid 
planes in desirable locations, and data must be read from files at each cross 
plane which could effect the thru- put time on the computer. 

A third possiblity is to switch to computing in a transformed compu- 
tational space. This makes the marching algorithm straightforward but 
forfeits one of the major advantages of GIM — completely arbitrary 
geometries. 

Approach 2 was selected as the best compromise and also provides the 
ultimate capability of elliptic -parabolic switching discussed earlier. 

The classic algorithms for treating the parabolic pressure field 
were deemed unsatisfactory. The following idea evolved from this 
research. The approach is termed "Quasi-Parabolic” and arose from 
the requirement of eliminating the ambiguities and singularities of 
existing methods. 
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The basic idea is to combine the classical parabolic marching approach 
with a *' quasi time’’ relaxation. The parabolic -march procedure greatly re- 
duces the amount of computer storage compared to a fully elliptic field. The 
time relaxation form of the equations eliminates the "decode” ambiguity asso- 
ciated with the parabolic pressure problem and allows velocity boundary con- 
ditions at solid walls to be treated. The equations used in the QP method are 
the time -averaged full Navier -Stokes, but with all second order terms dropped 
in a qua si -marching coordinate. Another way to view the QP equations is to 
take the parabolized Navier-Stokes and add back "psuedo time" derivatives. 

The QP solution procedure, as any parabolic marcher, thus allows no down- 
stream diffusion effects or pressure wave feedback through a subsonic flow. 

The solution is assumed known at upstream data planes, 1, 2, , . ,N-1, and the 
solution is sought at plane N with knowledge of plane N + 1. "Psuedo Time" 
relaxation, is used to obtain the solution at plane N in terms of the (converged) 
solution at a number of upstream data planes. Backward differences of some 
type, (second order) must be used to prohibit downstream feedback. So the 
QP algorithm is not a classical space marching scheme, and is also not a 
time -dependent elliptic method. It is somewhat of a hybrid technique which 
combines the better features of two approaches and eliminates the bad ones. 

The GIM/STAR elliptic code will converge a case in 500 to 1000 steps 
if the initial guess is chosen reasonably close to the answer. Also, GIM/STAR 
is relatively cheap to run, if the problem size is small enough to fit into memory 
and not require large page faults. The QP algorithm relieves both of these 
difficulties to some extent. By storing only a small number of data planes 
(and not the entire elliptic field) the large page fault problem is gone. The 
QP marching procedure can also assign a reasonable guess to the data 
plane since it knows the upstream converged solution, i.e., guess it is equal 
to the N-1®^ plane or extrapolated in some way. This should allow the time 
relaxation to converge very rapidly. If an implicit tlme-relaxer is used 
with the QP algorithm (with steps many times the CFL),the relaxation should 
go even faster. 
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The QP method allows an exact treatment of the Parabolic pressure 
field. No ambiguity exist in the QP decode at Mach = 1 (since it is '’quasi- 
time'* dependent) and no- slip walls can be treated exactly, i.e., boundary 
layers. The QP algorithm eliminates many of the bad features of pure 
parabolic methods. 

One obvious disadvantage of the QP approach is the planewise iteration 
(time relaxation) which must be done. This can be time consuming on the 
machine, and a good criterion for convergence must be used to avoid error 
propagation downstream. Spalding's method suffers from this same plane- 
wise iteration to correct the pressure as well; other linearization schemes 
such as Roberts and Forester (Ref. 5) which use the conservative equations 
also suffer. Planewise iteration is not uncommon in most parabolic methods, 
thus the QP scheme is no better or worse in this respect. A linearized block 
implicit scheme, as discussed in Section 3, appears to be very attractive for 
performing the quasi-time relaxation. 

Figure 4-1 shows the QP form of the three-dimensional Navier -Stokes 
equations in Cartesian coordinates. Note that these are the classical para- 
bolized form plus a psuedo-time derivative. Included are global mass con- 
servation, three components of momentum conservation, total energy and 
an equation for conservation of individual species in a binary mixture. 

Figure 4-2 is a typical computation molecule for a QP type marching. 

Assume that all flow variables are known at planes 1, 2, , .K and the solution 
is sought at plane K+ 1. If backward differences in x are used, (first order, 
second order, etc.), then the scheme of Fig. 4-3 allows no downstream 
feedback, and allows plane K + 1 to be uniquely determined from upstream 
information, i.e., qua si -parabolic ally. 

Now consider the ultimate, not immediate, implications of such an 
algorithm. A flow field could be marched out qua si-par abolically from an 
initial data plane 1 to plane K, where K is set a priori by the users. An im- 
bedded elliptic region is encountered between planes K and K+M. The 
number of planes that the QP algorithm can treat on any given sweep is not 
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Fig. 4-3 - The Quasl-Parabollc Scheme Finite Difference Equations. 



restricted to one . Simply specify single plane marching up to plane K, 
switch to the elliptic operator on the next M planes, and then para- 
bolically march from the (K+M) plane to the final N plane. The switching 
can then be done "automatically," but the user must still determine the loca- 
tion to perform the switching. Eventually, perhaps, an algorithm could be 
written to detect the onset of a separation bubble, flow reversal, or other 
elliptic phenomena. This is not being considered at this time, but only the 
fact that the capability is within the framework of the QP algorithm. 


Advantages of the QP Algorithm Outlined 


• There is no special treatment required of the parabolic 
pressure field. It is handled exactly except for the usual 
assumption of no downstream feedback. 

• No ambiguity exists in the decode procedure at Mach = 1. 
Thus mixed flows can be treated with no a priori knowledge 
of the relative velocity magnitudes, 

• Solid wall boundaries can be handled in the QP method 
with no -slip values. Regular parabolic procedures 
must avoid these type boundaries. 

• Inclusion of more than one upstream plane will allow 
second order accuracy to be maintained in the quasi- 
marching coordinate, 

• The QP scheme can accommodate either explicit or 
implicit "time" relaxation finite-differences. 

• Within the basic framework of the QP scheme, an 
elliptic region could be treated before, during or 
after a marching integration simply by including 
k-data planes (instead of 2) during the relaxation. 

• The QP algorithm requires very little addition stor- 
age over a classical parabolic method; and requires 
many times less storage than a fully elliptic treat- 
ment. Thus on STAR, the GIM/QP code could 
march out very large flow fields with no large page 
faults . 

• By dropping the cross -plane viscous terms, the QP 
procedure becomes Qua si -hyperbolic with no further 
coding changes. Thus one algorithm can accommodate 
either parabolic or hyperbolic flows. 
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4.4 RESULTS OF COMPUTATION 


The QP code has been essentially completed and a number of test - 
cases exercised. Three test problems are shown in this section for illu- 
stration of the Quasi-Parabolic code. These cases are: 

1. Flow in a three-dimensional duct with an expansion- 
recompression and interaction of two shock waves. 

2. Flow over a 10 degree planar wedge. 

3. Two-dimensional viscous flow resulting from interaction 
of a nozzle exhaust with a supersonic freestream. 

Other cases are currently in progress, including a boundary layer calculation, 
containing subsonic and supersonic flow. 

The first problem shown is depicted in Fig. 4-4. The 1x1 square nozzle 
expands via a trignometric variation to 10 units and then has a constant 2x2 
cross section. The supersonic Mach number expands up to the 10 unit plane 
then, due to recompression, shock sheets form at the top and outer side wall. 
The two shocks intersect as depicted in the figure. The QP code was used 
essentially in a hyperbolic mode with free slip inviscid solid walls and con- 
tained approximately 24,000 grid points. The purpose of the case is to deter- 
mine the ability of the QP marcher to handle rapid expansions and strong 
compressions (shock capturing). 

The solution for Mach number at the lower wall corner is shown on 
Fig. 4-4. A comparison is attempted here with a forward -marching hyper- 
bolic code of the classical variety (a MacCormack code). The GIM/QP solu- 
tion shows a strong shock wave while the other marcher would not solve for 
the large gradients at all. Figure 4-5 shows additional profiles for this case. 
The pressure ratio (local to inlet) is shown for both the upper and lower wall 
corners. Comparison is made to a published solution (Ref. 11). Excellent 
agreement is seen for the smooth upper wall profiles and for the expansion 
portion of the lower wall corner. At the axial location where the shock inter- 
section occurs, the two solutions differ considerably. The GIM/QP code, using 
a first order finite difference scheme agreed very well with the ATL results. 
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Fig. 4-5 - Three-Dimensional Square Duct GIM Hyperbolic Computation 
(Axial Pressure Distribution). 


However, the second order QP algorithm produced the curve shown in Fig. 4-5, 
i.e., a larger pressure rise. As a check on the accuracy of the QP shock wave, 
15,000 grid points were placed between l6 and 20 units. Very similar results 
were obtained as with the coarser mesh (11 x 11 x 81). It is thus felt that the 
GIM QP code is calculating the correct pressure rise across the shock. 

In order to test the shock-capture capabilities of the QP finite differ- 
ence scheme, an oblique shock on a 10-degree two-dimensional wedge was 
computed. Two example cases were run with incident Mach numbers of 1.8 
and 2.4. The same 60 x 51 node grid was used for both calculations. Each 
case required about 26 seconds to converge. The results are shown in Fig. 

4-6 as the pressure ratio through the shock as a function of vertical position 
and pressure rise from the NACA 1135 shock tables. The shock was char- 
acteristically smeared over five grid points. The excellent agreement indi- 
cates a good shock-capturing capability with the QP second order backward 
difference scheme. 

Case three consists of a parabolic, viscous flow in the configuration of 
Fig. 4-7. A nozzle with high pressure exhausts into a lower pressure, hyper- 
sonic freestream flow. This case was solved with the GIM elliptic code with 
940 nodes and reported in Ref. 2. The QP algorithm gives virtually identical 
results as given by the full Navier -Stokes code. The grid used and the steady 
state Mach and pressure contours are shown on Fig, 4-7. Comparison of this 
solution with the reported values of Ref. 2 and with the inviscid SEAGULLi code 
of Ref. 12 are shown in Figs. 4-8 and 4-9. The SEAGULL is an inviscid, slip- 
line, shock fitting, forward marching code. Figure 4-8 shows vertical pres- 
sure distributions at three axial stations in the shear region, and Fig. 4-9 
gives the corresponding Mach number plots. As seen by the comparisons, 
the GIM marching algorithm does indeed work as expected and gives quanti- 
tatively the same answers as the other codes. Application to a boundary layer 
problem is currently under way. 
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Fig. 4-6b - Wedge Shock Case to Veri 
(M = 2.4, 6 = 33 deg). 
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Fig. 4-7 - Quasi-ParabolLc Code Applied to Shear Flow 
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Fig. 4-8 - Quasi-Parabolic Shear Flow Computation (Vertical Pressure Distributions 
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Appendix A 


A.l THE GIM/STAR SE-2 CODE 

The Blue Book (Ref. 1) describes the version of the GIM/STAR code 
designated SE- 1 (STAR Elliptic Version 1). This reference manual contains 
input guides and user information for the code. Since the publication of this 
Blue Book there have been a number of changes to the code which have not 
been documented. Some of these changes were necessary to allow large 
problems to be run with a minimal number of large page faults while others 
were made to reduce the possibility of wasting computer time generating 
MATRIX analogs on a bad grid. 

The input changes are not extensive but the user should use this 
Appendix in conjunction with the Blue Book when running a GIM/STAR 
problem. The following subsections describe the changes for the program 
modules and file usage. 

A. 2 GEOM MODULE 

Module 1 of the GIM SE-1 deck was titled GEOMAT as it contained both 
the geometry and grid generation and the matrix coefficient assembly. The 
SE-2 version has the two operations broken out into separate modules. The 
first module of SE-2 is titled GEOM, as it now only performs the geometric 
description and grid generation (see Fig. A-1). 

The user should be aware of the differences in this module between 
versions 1 and 2: 

• Input cards l6 and 17 (in the Blue Book) are no longer used — just 

omit cards 16 and 17. 
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Card Type 


Parameter Liist/Format 


1 HEADER(I), I = 1, 72 

(12A6) 

2 NZONES, IDIM, ISTEP, IMATRX, IMATE 

(615) 

3 IWRITE, LWRITE, NWRITE 

(315) 

4 KC(I), 1=1,6 

(6A5) 

5 NSECTS 

(15) 

6 MAPE(I), 1=1,12 

(1215) 

7 MAPS(I), I = 1,6 

(615) 

8 (IBWL(I), I = 1,6), ITRAIN 

(715) 

9 (NNOD(I), 1=1, 3), (ISTRCH(I), 1=1,3) 

(615) 

10 DIVPI(I), I = 1, 3 

(3E10.4) 

11 [AETA(J,I), 1=1, NNOD(J)] , J = 1, IDIM 

(8E10.4) 

12 [(AC(I, K, J), I = 1,8), J = 1,4 or 12], K = 1,5 

(8E10.4) 

13 [AS(I, J), I = 1,8], J = 1,6 

(8E10.4) 

14 (PT(I, J), I = 1,5), J = 1,4 or 12 

(8E10.4) 

15 [(PMAX(I,K, J), I = 1,5), ETAMAX(K, J), K = 1,4] , 
J = 1,4 or 12 

(6E10.4) 


Flg.A-l - Input Guide for GEOM Module (SE-2 Code), 
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• File 17 is not output from the geometry module, rather File 18 is 
now to be saved. This new File 18 is to be subsequently input to 
the new MATRIX module. 

• Card Type 4 has been changed. The values of a are no longer 
input, but rather a set of flags to retrieve the correct a’ s are 
now used. The parameter KC is set to alphabetic characters 
F, B, etc., for forward or backward differences, KC(1) is for 
X step 1, KC(2) for y step 1, etc,, through KC(6) which is z 
step 2, The format is 6A5. 

A. 3 MATRIX MODULE 

The new MATRIX module of the SE-2 code now performs the analog 
coefficient calculation and file creation. This module should be executed 
following GEOM and before INTEG. The File 18 which was output from GEOM 
is now input to the MATRIX module. File 17 needed by INTEG is to be saved 
from MATRIX. Figure A-2 gives the storage requirements for MATRIX. 


The input to the MATRIX module consists of the Cards 16 and 17 which 
were omitted from the GEOM module, plus one new card. Each of these cards 
is now described: 


Card 

Parameters 

Format 

1 

NDX, NDY, NDZ, ISNOPT 

(415) 

2 

KC(I), 1= 1,6 

(6A5) 

3 

Nl, IC, NT 

(315) 


Card Type 1 Format (415) 

Same as Card Type l6 (GEOMSE-1) p. 4-27 


nodal decrement in the rj coordinate system 


NDX 

NDY nodal decrement in the rj system 

NDZ nodal decrement in the system 

ISNOPT special node treatment flag 

If ISNOPT = 1 the MATRIX module will calculate the 
number of special node terms placed 
on File 17 for input to INTEG 
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Matrix 

2D 

3D 

/A COM/ 

50*NN 

196=?‘NN 

/pcom/ 

4^NN 


(Q3MAP/ 

24*NN + IS^SPEC 

48*NN + 18*NSPEC 


+ 65542 

+ 65542 

/IRFBC/ 

6*NSPEC 

6 ❖NSPEC 

/JCFBC/ 

6'1«NSPEC 

6*NSPEC 

/PAFBC/ 

6 -’‘NSPEC 

6 ❖NS PEC 

where 



NN 

= total number of nodes 


NSPEC 

= number of special node terms allowed 
for in DIMENSION statements (DYNMAT 
input) 

The common block sizes may be calculated for each problem 
size to determine the ideal grouping on the LOAD card. If in 
doubt assign each block to a new large page boundary as below. 
Do not use GRLPALL, but use 

GRLP=^ACOM, GRLP=^^PCOM, GRLP=-*Q3MAP, GRLP=IRFBC, 

etc. 




Fig.A-2 - Module Common Block Sizes. 
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If ISNOPT = 0 


the entire array of special node terms 
will be placed on File 17. The size of 
the array is determined by DYNMAT 
input (NSPEC parameter). 


Card Type 2 Format (6A5) Analog Choice Card 

This card consists of a sequence of six characters (F or B) 
identifying the difference direction (forward or backward) 
for X, Y, and Z Step 1 and X, Y, Z Step 2, respectively. 


Examples: 


F F F B B B 


forward, forward, forward, backward, backward, backward 
for three-dimensional problems and 

F B _ B F _ 

forward, backward, backward, forward for two-dimensions. 


This is a new card for version SE-2 and is identical to GEOM card type 4. 


Card Type 3 Format (315) 

Nodal analog print control card. 

N1 first node of a print sequence 

IC print increment 

NT total number of nodes to print for this sequence. 
(See page 4-29 SE-1 manual for complete description.) 

Any number of cards of this type may be input. 

Place a - 1 in Columns 4 and 5 on last card to terminate. 


Dynamic Dimension for MATRIX 


The new MATRIX module has its own dynamic dimensioning sequence. 
The same deck DYNMAT is used for MATRIX and GEOM, but the input of a 
third parameter is optional in MATRIX. 
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deck. 


The DYNMAT deck is to be executed before running the MATRIX 


The input consists of one card; 


Card 


Parameters 


Format 


1 


NN, IDIM, NSPEC 


(315) 


The definition of these input variables are the following; 


Card 1 Format (315) 

NN number of nodes 

IDIM dimensionality (2 or 3) 

NSPEC number of special node terms to allow for in 

DIMENSION statements. If left blank, or zero, 
the arrays will be dimensioned to NN. The 
actual number of special node terms will be 
calculated and printed out in MATRIX. This 
value is then input to INTEG. (Not used by 
GEOM.) 


A. 4 INTEG MODULE (SE-2) 


This module has remained virtually intact from a user standpoint. 
Three additional options have been added since the Blue Book was issued. 
These are; 

• Capability to compute a CFL time step automatically over 
a multi -zoned grid. 

• Treatment of downstream subsonic boundary conditions using 

a mass balance condition. This option was added under another 
NASA contract and is documented here for completeness. 

• Input of a set of flags denoting the finite -difference direction. 

This aids in a more complete set of difference options and allows 
for full vectorization of all schemes. 
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Figure A-3 is a summary Input Guide for the INTEG SE-2 module. 

Note that the new input cards are designated 2a, 2b, 2c, 3a and 6a. All except 
6a are optional and existing data decks will still work as they did for version 
SE-1. Card 6a must now be input in version SE-2. Card 2 has two additional 
inputs, IDS, IBOUND which control the optional input of 2a, 2b, and 2c. Zero 
values for the parameters on Card 2 signifies omission of the remaining Cards 
Type 2a, 2b and 2c, Figure A-3 is a description of the available options and 
each parameter that is to be input. 

Figure A -5 describes each parameter that is input on the optional Cards 
2a, 2b and 2c. An example of the use of the subsonic boundary condition option 
is shown in the sample grid of Fig.A-6, 

The time step calculation option is controlled via the value of KZONES 
read as the last data on Card 3. If this is omitted (=0), then one zone is 
assumed. If KZONES > 0, then this signals the code that a multiple zone 
problem is being run. In this case, the value of KZONES should be equal to 
the number of zones used in the geometry module. If KZONES = 0, then Card 
Types 3a are not used, but any value of KZONES > 0, requires the input of 
come Cards 3a. The number of Cards 3a to be input is equal to KZONES- 1. 
The time step information for zone number 1 is input on Card 3 itself. 

Figure A -7 describes the input of this time step information. 

Card 6a is simply the KC values used in GEOM and MATRIX, i.e., 

, F F F B B B 

in format 6A5, This card must agree with the previous module’s usage. 

One additional Fig. A -8 is included in this subsection. This chart 
shows formulas for determining the COMMON block sizes for the INTEG 
module. These values are needed for large problems to set up the LOAD 
card as described in the Blue Book (Ref. 1). 
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Card Type 
1 

2 


2a 


2b 


2c 


3 


3a 


4 

5 

6 


6a 


7 

8 

9 

10 
11 
12 


Parameter Llat/Format 

ICASE, IITITLE(I), 1=1, 78) 

(12, 78A1) 

IDIM, METHOD, ITMAX, IPRNT, ITSAVE, ISTART, 
lOTYPE, lUNITS, ITSTRT, IVISC, IDIST, ISPEC, 
IDS, IBOUND 

(1415) 

INFOUT, IJUMPO, JJUMPO, NIOUT, NJOUT, 
ICALC, AMFLW 

(615, ElO.O) 

INFINL, IJUMPI, JJUMPI, NUN, NJIN, ICALC, 
OUTMFL 

(6I5, ElO.O) 

INFINL, IJUMPI, JJUMPI, NUN, NJIN, INFOUT, 
IJUMPO, JJUMPO, NIOUT, NJOUT, ICALC 

(1115) 

NN, NNX, NDX, NNY, NDY, NNZ, NDZ, NPM, 

K ZONES 

(915) 

KST, KNX, KDX, KNY, KDY, KNZ, KDZ 
(715) 

DTIME, DTFAC, INCDT 
(2E10.0, 15) 

REALMU, REALK, GAMSl, GAMS2, WMl, WM2, 
DK, RK 

(8E10.0) 

EMU, ELAM, ERHO, ESPEC 
(4E10.0) 

KC(I), 1= 1, 6 
(6A5) 

NNPM(I), NCPM(I), (NNCPM(I, J), J = 1,5), 
ANGPM(I); I = 1,NPM 

(715, ElO.O) 

(NCT(I,J,K), PXPM(I,J,K), PYPM(I, J,K), 

K = 1,4); J = 1, NCPM(I); I = 1, NPM 

(15, 2E10.0) 

RHOZ, PZ, ASTAR, NINC, A, B 
(3E10.0, 15, 2E10.0) 

NJ, INC, NTOT, ITAN, ITYPE 
(515) 

RI, UI, VI, WI, PI, CSI 
(6E10.0) 

Nl, IC, NT 
(315) 
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Fig. A-3 - Input Guide for INTEG Module (SE-2 Code) 



I 


Card 

Col. 

Format 

Variable 


Description 

Type 2 

1-5 

15 

IDIM 

See Blue Book 


6- 10 



METHOD 

1 



11-15 



ITMAX 

! 



16-20 



IPRNT 




21-25 



ITSAVE 

j 



26-30 



ISTART 




31-35 



lOTYPE 




36-40 



lUNITS 




41-45 



ITSTRT 




46-50 



IV ISC 




51-55 



ID 1ST 

i 



56-60 



IS PEC 


r 


61-65 

1 


IDS 

Boundary Condition Flag 






= 0, one-sided differences 






(supersonic) 






= 1, mass balance technique 






(subsonic) 


65-70 

15 

IBOUND 

Note; 

If IDS. Eq. 1 IBOUND 


should be set to either -1, 0, or 1. 

If IDS. Eq. 0, IBOUND Is left blank 

= -1, input inlet mass flow and 
calculate exit mass flow 

= 0, input exit mass flow and 

calculate inlet mass flow 

= 1, calculate both inlet and exit 

mass flow 

Note; If IDS. Eq. 1 card types 2a, 2b and 2c must follow type 2 card. The use 
of types 2a, 2b and 2c depends on the value of IBOUND. 

If IBOUND = -1, use Type 2a 

If IBOUND = 0, use Type 2b 

If IBOUND = 1, use Type 2c 


Fig. A -4 - Definition of parameters for Card 2. 
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Card 


Type 2a 


Type 2b 


Type 2c 


Col. 

Format 

Variable 

1-5 

15 

INFOUT 

6-10 

15 

IJUMPO 

11-15 

15 

JJUMPO 

16-20 

15 

NIOUT 

21-25 

15 

NJOUT 

26-30 

15 

ICALC 


31-40 

ElO.O 

AMFLW 

1-5 

15 

INFINL 

6-10 

15 

IJUMPI 

11-15 

15 

JJUMPI 

16-20 

15 

NUN 

21-25 

15 

NJIN 

26-30 

15 

ICALC 

31-40 

ElO.O 

OUTMFL 

1-5 

15 

INFINL 

6-10 

15 

IJUMPI 

11-15 

15 

JJUMPI 

16-20 

15 

NUN 

21-25 

15 

NJIM 

26-30 

15 

INFOUT 

31-35 

15 

IJUMPO 

36-40 

15 

jjUMPO 

41-45 

15 

NIOUT 

46-50 

15 

NJOUT 

51-55 

15 

ICALC 


Description 

Starting node on exit plane 
th 

Nodal increment in i direction on exit 
plane 

Nodal increment in direction on exit 
plane 

Number of elements in i^^ direction on 
exit plane 

Number of elements in j direction on 
exit plane 

Velocity update flag 

= 1, update inlet velocities 

= 2, update exit velocities 

Inlet mass flow rate (input by user) 

Starting node on inlet plane 

th. 

Nodal increment in the i*^ direction on 
inlet plane 

Nodal increment in the direction on 
inlet plane 

th 

Number of elements in i direction on 
inlet plane 

th 

Number of elements in j direction on 
inlet plane 

Velocity update flag (see Card Type 2a) 

Exit mass flow rate (input by user) 

Starting node on inlet plane 

th 

Nodal increment in the i direction on 
inlet plane 

th 

Nodal increment in the j direction on 
inlet plane 

Number of elements in i*"^ direction on 
inlet plane 

th 

Number of elements in j direction on 
inlet plane 

Starting node on exit plane 

th 

Nodal increment in the i direction on 
exit plane 

th 

Nodal increment in the j direction on 
exit plane 

Number of elements in i^ direction on 
exit plane 

Number of elements in direction on 
exit plane 

Velocity update flag (see Card Type 2a) 


Fig. A -5 - Description of Input Parameters for Optional Card Types 2a, 2b, 2c (Subsonic 
Boundary Conditions). 
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6 18 30 42 54 66 

68 
70 
72 
74 

J=1 J=2 J=3 J=4 J=5 


1= 1 

8 

20 

32 

44 

56 

1 = 2 

10 

22 

34 

46 

58 

1=3 

12 

24 

36 

48 

60 

1 = 4 

14 

26 

28 

50 

62 


Example; 

IBOUND = 0 INFINL = 6, IJUMPI = 2, JJUMPI = 12, NUN = 4, NJIN = 5 


Fig. A-6 - Example of Subsonic Boundary Condition Usage. 
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Card 


Description 


3 


Col. 

Format 

Parameter 

5 

15 

NN 

10 

15 

NNX 

15 

15 

NDX 

20 

15 

NNY 

25 

15 

NDY 

30 

15 

NNZ 

35 

15 

NDZ 

40 

15 

NPM 

45 

15 

KZONES 




See Blue Book 


The number of zones that was used to construct the grid. This is used to allow 
a CFL time step to be computed over multiple zones. Set to 1 for a single zone 
problem. 


Card 

Col. 

Format 

parameter 

Description 

3a 

5 

15 

KST 

Starting node number 
of this zone. 


10 

15 

KNX 

Number of nodes in rj ^ 
direction for this zone^ 


15 

15 

KDX 

Nodal decrement in 
direction for this zone 


20 

15 

KNY 

Number of nodes in 
direction for this zone 


25 

15 

KDY 

Nodal decrement in TJ ^ 
direction for this zone 


30 

15 

KNZ 

Number of nodes in 
direction for this zone. 
Set to 1 for 2“D flow. 


35 

15 

KDZ 

Nodal decrement in 7J^ 


direction for this zone. 
Set to 1 for 2-D flow. 


Note; Input Card Type 3a for each multiple zone to be used in computing a CFL 

time step. The number of Cards 3a is equal to KZONES-1, where KZONES 
is input on Card 3. 


Fig. A -7 - Description of Parameters for Optional Card Type 3a Input. 
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Common Block 
Names 

Axisymmetric 

2 

-D 

3- 

-D 

1 Gas 

2 Gases 

1 Gas 

2 Gases 

1 Gas 

2 Gases 

/prim/ 

5*NN+1 

6^«NN+1 

5*NN+1 

6*NN+1 

6'i^NN 

7^'NN 

/TAU/ 

9’^NN+3 

9-^NN+3 

9*NN+3 

9-NN+3 

12*NN 

12*NN 

/TMVEC/ 

2«NN 

2*NN 

2»NN 

2*NN 

2^i^NN 

2^NN 

/VPROP/ 

2«NN 

4^^NN 

2^NN 

4*NN 

2^^N 

4*NN 

/vbuf/ 

8^mN 

10*NN 

8«NN 

10*NN 

lO^l^NN 

12'^NN 

/BOUND/ 

5^NB + 1 

5^NB+1 

5»mB+l 

5*NB+1 

5^'NB+l 

5^NB+1 

/ebuf/ 

8 'i'NN+4 

10=J'NN+5 

8 >J«NN+4 

10*NN+5 

15^^NN 

18=J'NN 

/xbuf/ 

7*NN+4 

7*NN+4 

7 ^NN +4 

7 *NN+4 

10*NN+1 

10«NN+1 

/step/ 

3^NN+10 

3*NN+10 

3*NN+10 

3^^NN+10 

4 *NN+9 

4*NN+9 

/axsym/ 

8*NN 

9*:'NN 

8 

9 

9 

10 

/Q3MAP/ 

24^NN+ 18‘^^NSPEC 
+ 6+ COM!P 

24 *NN+ 1 8 *NS PE C 
+ 6+ COMP 

48>:^NN+18*NSPEC 
+ 6+COMP 


NN 

NB 

NSPEC 

COMP 


= total number of nodes. 

= number of boundary nodes. 

= number of special nodes. 

= amount of storage need to complete a large page. 


Fig. A-8 - INTEG Module Common Block Sizes for SE-2 Code. 
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A.5 GIM SE-2 FILE DESCRIPTIONS 


Following is a brief description of the files used in the STAR SE-2 
system. In all but very small problems setups, (a few hundred nodes), a 
REQUEST card must be used for each file. The form of the REQUEST 
card is as follows: 

REQUEST (FILEXX/NSPGS, T = P) 


where 


NSPGS = the number of small pages of disk space 
allocated to the file 

1 small page = 512 words 


Formulas for calculating NSPGS are now given for each file. 

FILE 16 GEOM 

Work file used in GEOM only (Binary) 

NSPGS ^ 15=^'NN/512 NN = no. of nodes 

FILE 17 MATRIX/INTEG (Binary) 

Nodal analog file created in MATRIX and used in INTEG. 


NLPGS = (16*NN + IS^NSP + 6)/65536 

rounded up to next whole integer 

NSPGS = NLPGS *128 + 1 


3D 


NLPGS = (48*NN + 18*NSP + 6)/65536 

rounded up to next whole integer 

NSPGS = NLPGS *128 + 1 


where 

NN = total number of nodes 
NSP = number of special node terms. 
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FILE 18 GEOM/MATRIX (Binary) 


File containing matrix assembly data. Created in GEOM and used in 
MATRIX. 

2D 

NSPGS = 50^-NN/512 + 1 
3D 

NSPGS = 196»NN/512 + 1 

FILE 20 GEOM/INTEG/GIMTEK (Formatted) 

Nodal geometry file created in GEOM and used in INTEG and GIMTEK 
2D 

NSPGS ~ 14*NN/512 
3D 

NSPGS ~ 20-^NN/5l2 

FILE 22 INTEG/GIMPLT (Formatted) 

Flowfield solution file created in INTEG and used both as a restart 
file and in GIMTEK. 

2D 

NSPGS ~ 11^5'NN/512 per record 

NSPGS ~ 14-'=NN/512 per record 

Multiply by the number of iteration increments saved. 
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Controllee File Sizes 


The size of the controller file is specified on the LOAD card in small 
pages. Formulas are given below for calculating the size required for a 
given problem. 


GEQM 

( 50WN/65536 + 3 2D 
NLPGS =< 

(196^J^NN/65536 + 3 3D 
rounded up to next whole integer 


NSPGS = NLPGS*128 This is the value that goes 

on the LOAD card. 


MATRIX 


No single formula exists to calculate the controllee file size for the 
MATRIX module. The procedure is to calculate the number of large pages 
(65536 words each) required for each GRLP parameter on the LOAD card, 
add these up, add 2 for other storage and multiply the result by 128. 

INTEG 


The same rule applies to the INTEG controllee file size as to the 
MATRIX module. Use the common block sizes to compute the number of 
large pages, add them up, add a couple and then multiply by 128 to get the 
controllee file size number. 

A. 6 PLOT MODULE (GIMTEK) 

The GIM SE-2 code plotting module is now titled GIMTEK. This re- 
flects the modifications which were made to the GIMPLT SE-1 module in order 
to use the Tektronix 4014 for graphic output. The user need not be aware of 
the internal program changes that were made. The input data are identical to 
the SE-1 version. Three items of significance to the user are now described: 
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KMAX^q = max 


5040 + 8*NN 

1 1 ^'NN 


NN = number of nodes 


The CM field length specified on the job card is calculated by 


CMjq = KMAX^q + 23000^q 


Notes: 1. This parameter must be set in the program and the array "A 

dimensioned to this value. 

2. CM must be converted to octal for specification on job card, 
and RFLj card 

3. Example 


NN = 2000 


KMAXjq = max 


5040 + 8*2000 


11*2000 


= max 


21040 


22000 


= 22000 


CM^q = 22000 + 23000 = 45000^^ 


use CM = 130000 (130000o = 45056 , ^ > 45000 ^ 


Fig. A -9 " GIMTEK Core Requirements 




• The formulas on page 6-25 of the Blue Book for computing core 
sizes for the plat module are no longer valid. Figure A -9 gives 
the revised formulas and an example calculation. 

• The plot save command was changed on the software system. The 
new save name is 


SAVPVF. 


• The routine that we use for obtaining GIMTEK plots from the 
Tektronix 4014 is 


PLIST. 

This allows enough options to select only those plots needed and 
also allows an unlimited time to examine a plot before proceeding. 


The input data for GIMTEK is the same as described in the Blue Book. 
Figure A -10 is a summary of the required input data presented here for 
completeness. The user is referred to the Blue Book for a definition of the 
parameters. 


A. 7 EXAMPLE RUNSTREAMS FOR THE SE-2 CODE 

The following pages show example runstreams that have been used for 
the SE-2 code on the STAR- 100 machine. These should aid the new user in 
setting up a deck for GIM SE-2: 

Fig. A- 11 - GEOM Module Only 

Fig. A- 12 - MATRIX Module Only 

Fig. A- 13 - GEOM/ MATRIX Combination Run 

Fig. A- 14 - INTEG Run Only 

Fig. A- 15 - GIMTEK Run 


Note: The "blanks” which show up on the card listings are 7-8-9 punch cards, 
i.e., end of record. 

The files for GIM SE-2 are cataloged under user number 838700C as GEOMB, 
MATRIXB, INTEGB and GIMTEKB. 
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Card Type 

Parameter List/Format 

1 

ITITLE(l), ITITLE(2) 
(2A40) 

2 

NX, ITERAD, ITRBLK, KDIM. ISP 
(515) 

3 

GAMMA, FACTOR, RK, PO, TO, RHOO 
(6E10.0) 

Specs. 

S-1 

NPLT, STITLE, IVIEW, ISYM, ITHETl, lAXISl, 
ITHET2, IAXIS2, IXTABL, VFAC 

(15, 5X, A20, 815, ElO.O) 

S-2 

NTYPE, JO, IJUMP, JJUMP, NI, NJ, IPRNT 
(715) 

Grid 

G-1 

'GRID', lOPT, ICSCLE, NSPECS, (ISPEC(I), 1=1, 
NS PECS) 

(A4, IX, 15, 26X, 215, 715) 

G-2 

(ISPEC(I), 1 = 8, NSPECS) if NSPECS > 7) 
(45X, 715) 

VVEC 

V-1 

'VVEC, lOPT, NITER, ICSCLE, NSPECS, (ISPEC(I), 
1=1, NSPECS) 

(A4, IX, 215, 20X, 215, 715) 

V-2 

(ISPEC(I), 1 = 8, NSPECS) (if NSPECS > 7) 
(45X, 715) 

I-l 

(ITER(I), 1=1, NITER) 
(16I5) 

Contours 

C~1 

ITYPE, lOPT, NITER, NC, ITABLE, INCR, ICSCLE, 
NSPECS, ISPEC(I), 1=1, NSPECS) 

(A4, IX, 515, 5X, 215, 715) 

C-2 

(ISPEC(I), 1=1, NSPECS) (if NSPECS > 7) 
(45X, 715) 

I-l 

(ITER (I), 1=1, NITER) 
(1615) 

L-1 

(CVAL(I), 1=1, NC) 


(8E10.0) 


Fig. A- 10 - GIMTEK Summary Input Guide. 



GEuDCT « CM60000 * T 1 OO • 

USE;h:,838700C. 

CHAkGE* 101857*LKC. 

GET( OLDPL=GEOM/UN=838700C) 

GET ( DYNmaT=DYNiviAT/UiM = 838700C > 
UPDATE ( F , C=TAPE8 ) 

DYNi'-iAT. 

T0 STAR( InPUT^TAPEs) 

*ID MODS 
254 1 3 

STUPE 838700 400SDS TESTDECK 3 
STRSIOE^TIOO. 

FOPThiAN( I = TAPE3<B = GE0Mb»O = LS) 
PEUUEST(fILE16/75, T=P ) 

REOUEST( file 18/974 *T = P ) 
request (F ILE20/ 1 00 » T =P ) 

LOAD ( GEOMB f CN^GEOiMGU 9 1408 » GRLPALL 
GEOi^IGO, 

TO AS ( Z =838700C* F ILE18 = B I »F ILE20 ) 
*** GEOMETRY DATA **■»■** 


Fig. A- 11 - GEOM Runstream 



I 


MATRIX»CM6000U*T100. 

USER f 63fcl700C. 

CHARGE* lO 1857*ERC, 

GET t Oi_UPL = MATRJ X/uN = 838700C ) 

GET( UYNi-iAT = DYf\MAT/UN=838700C ) 

UPDATE (F, C = TAPE8 > 

COPYSLiF( tapes* OUTPUT ) 

REWI ND< TAPES) 

DYNMAT • 

ATTACH(FILE18=FI LEibB ) 

TOb I ARC INPUT * TAPES >F I LE 18 = 8 I//.U) 

* I D NONE 

Sb41 3 1714 

STuNE 838700 400SDC3 TtSTDECK B 
STRbIDE*T 100. 

FOR 1 RAN( I =TAPe3*B = MATRB*0 = LB> 

REOUEST (FI EE 17/385* T = P ) 
load ( MATRB*CN=rvlATRGO * 1920 

iGREP=:->^AcQM* GREP = *0 3)’|AP ♦ GRLP = 'J^PCUi''i * ■>«■ 1 RF BC * JCFBC * *PAF oC ) 
/viATRGO. 

TOAb(Z=B38700C*FlLE17=BI ) 

121 11 1 i 

F F F d B B 

1 1 20 

- 1 


Fig. A- 12 - MATRIX Runstream. 
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GEC;|'''iaT * CM6000U » T 1 OO • 

USEH^*63fc5700C. 

CHAh-GE* 101857, LRC, 

GET ( ULOPl = GEG-^i/UN = 83870OC ) 

GET t DYNiv'iaT = DYNfVlAT/UiM = 838700C ) 

UPDATE ( h, C=TAHE8 ) 

DYNi'/IAT • 

C0PYCF(TAPE3,GEC vic) 

PEWlNDt G eOMC) 

RETURN( OLDPL ) 

RETURN( TaPE3 ) 
return ( Tapes ) 

GET ( OLDPL = .’^ATrxI X/U)V| = 838700C ) 

UPDATE (F,C=TAPE8 ) 

dynmat* 

COPYCF( TaP£3*MATC ) 

REW I ND ( fv.ATC ) 

RETURN ( ULDPL ) 

RETURN( TAPE3) 

RETURN( TAPES) 

TOST AR ( I NPuT I GEUNc , .^lA f C ) 

* I D MODS 

254 1 3 

* I d none 

254 1 31714 

STU^<E 838700 400SDS TtSTDECK 8 
STRSIDE,T100« 

REOoEST (FI LE16/75, T = P } 

REUULST(F ILE 17/388 » T =P ) 

RE0^EST(F I LEI 8/9 74 » T =P ) 

REOUEST (F I LE20/1 OQ 9 I'=P ) 

FORTRAN ( I =GEOKlC, 3 = G£0MB » 0 = LB ) 

LOAD ( GEOMB, CN^GEOmGu , 1408 , GRLPALL= ) 

GEOjv.GO. 

F0RTRAN( I =MATC, B = ;vlATRt3 >0 = LB ) 
load ( MATRB , CN = ,ViATRGU , 1920 

, Gki-P = * AGON , G>^LP = 3i''lAP * GRLP = *PC0i''l I RFbC ,*JCFtjC» '^RAPbC ) 

i''lATRGO • 

TO AS ( Z =838700C, F I Lt 1 8 = 8 i’, F 1 LE 1 7 = B I ,F I LE20 ) 

***:«■ Geometry data 

MATRIX DATA ;t*^(-^f*** 


Fig. A- 13 - GEOM/MATRIX Runstream 



1 NT EGA » CM6000'J « T200 • 

USEt^f 012839C. 

CHA^^iGE^ lOH I 10»EKC. 

ATTACH (F I LEI 7 = F I LeI 7A ) 

ATTACH< F ILE20=F I Le20A J 

GET ( OLDFL= INTEG/UN=83a700C) 

get t DYNDIIv. = DYNDI H/UiM = 83ti700C ) 

UPDATE ( F, C=TAPE 8 ) 

DYND IM, 

COPYCFC TAPE3f INTEgX ) 

REW I ND < I NTEGX ) 

RETURN( OLDPL ) 

RETURN ( TAPE3 ) 
return ( Tapes > 

TUbTAR ( 1 NPUT ♦ i NTEgX >F 1 LE20 » K I LE 1 7 = Li 1 //U ) 

I U L N T MODS 

122 S 1 0 I7b9 

stoke 012839 AOO^DS TlcjTDECK B 
STROIDE»T200. 

FOk TRAN( I = I NTEGX , 3 = I f £GB/ 1 00 » 0=LB ) 

LOAD( I NTEGLi* CN= I NtEGO > 2000 
• GRi— P = •>''PR I |V' « -X- T ixlV EC > VF i\OP * T AU 
, *Ei-JUF » *UBUF ♦ abound > "'^'XyUF ♦^STEP 
» GkuL = '^'>-J3I'^1AP ) 

INTEGO, 

TO AS ( Z=Ol 2839C 1 F IlL22 ) 

***** INTEG data **^****** 


Fig. A -14 - INTEG Runstream. 


101 



GIMTEK*CM 1 20000* T 40 a . 

USER*012839C. 
charge* 1o21 10*LRC, 

GET ( OLDPl = GIMTEK/uiM = A92429C ) 

UPDATE ( F ) 

FTNt I =CO!VlP IEE’E=0 ) 

ATTACH( TAPE20=FILe2JA ) 

ATTACH ( TAP£22=F I LE22E ) 

PFL C 120000 ) 

ATTACH (i_ I bPTEK*LHCG05h /CTN^L IBPAPY) 

LDsE-T (i_lB = L I BKTEK/LRCGOSF , preset =:NG1 IMF ) 
LGO. 

SAVE ( SAVPVF=SCRJET ) 

*ID KORCHG 
^'<■1 GIiMPlT.744 

ISeT= lSET+1 
*D GIMPlT*9 

COMMON A( 14840) 

GIMPlT.15 
KMAX= 14S40 

^-Jt*** GIMTEK data -;f 

STOP 


Fig. A - 15 


GIMTEK Runstream 
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Appendix B 


B.l INTRODUCTION 

Algorithms are developed for the solution of the three-dimensional 
compressible Navier -Stokes equations in conservation form. This work 
represents an extension of the methodologies outlined by Beam and Warming 
(Ref. B-1) and Spradley (unpublished information) and familiarity with the 
cited literature is assumed. A time -dependent algorithm for the unsteady 
equations is developed first and then a spatial marching scheme for the 
three-dimensional parabolized steady equations is obtained. Algorithms 
for one- or two-dimensional problems are easily obtained by simply de- 
leting appropriate terms from the equations. 


105 


B.2 THREE-DIMENSIONAL UNSTEADY ALGORITHM 

The three-dimensional compressible conservation equations can be 
written in conservative form 

at ’*■ 8 x ay az ■ 

= + Vi2<U.Uy) + Vj3(U,U^)] 

+ ^["^21 (U. U^) + V 23 (U. Uy) + V 23 (U, U J] 

+ + V32(U.Uy) + V 33 (U,U^)] (B.l) 

where U is the vector of conserved variables and E, F, G and V.. are flux 
vectors. 


A generalized single-step temporal finite difference scheme for ad- 
vancing the solution of Eq. (B.l), is the following. 


AU 


n 


9 At aAu^ au” 

1+1 at 1+^ at i+e 


AU 


n-1 


(B.2) 


where U^ = U(nAt) and AU^ = U^"^^ - u'^. (Terms of order At^ and At^ have 
been neglected, for simplicity.) 
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au 


If Eq, (B.l) is solved for and inserted in Eq. (B.2), the resulting ex- 
pression for AU”" is 


AU 


n : 9 At 

" 1+e 


^ (-AE" + Av“j + Av"2 + Av“3) 


+ (-af“ + AV21 + AV22 + AV23) 


8y 


+ h + ^^32 + ^^ 33 > 


At 

1 + ^ 


3x 


(-E + Vji + Vj2 + Vjj) 


n 


+ ^ (-F+ ^21 + ''22+^23)” 


+ ^ ^ 31 + ^^ 32 + ^ 33 >" 


+ -=1^ AU*'-' (B.3) 

where AE^ = E^^^ - E , etc. 

Note that AE^, AF^, AG^ and AV^^ are nonlinear functions of the con- 
served variables U. A linear equation with the same temporal accuracy as 
Eq. (B.3) can be obtained by expanding AE^, AF , AG and in a Taylor 

series, thus 

I- . E” * 0 AU" 

or 

Ae"^ = AU’^ (B.4a) 
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Similarly 


AF 


n 


AG 


n 


(0 

(B,4b) 

(i)" 

(B.4c) 


where A, B, and C are the linearization Jacobians . 

Strictly speaking, the same procedure should be applied to the viscous 
terms. However, as pointed out in Ref. B.l, treating the spatial cros s -derivative 
terms, i.e., AV-. (i/j), in this manner would lead to considerable difficulties 
in constructing an efficient spatially factored algorithm. Therefore, spatial 
cross-derivative terms will be evaluated explicitly (without loss of accuracy, 

Ref. 1), i.e., 


AV" = Avl^."^ (i;^3) (B.5) 

while the linearization is applied to the (k = 1, Z, 3). Remembering that 

^^kk = f(U,U ), 


AV 


n 

kk 


av 


n 


kk 


BU 


. n 
AU + 


dv 

Wu 


n 


kk 


AU 


X, 


n 

X, 


= Pfck + Rkk 


(B.6) 


Application of the product differentiation rule shows that 


R 


n 

kk 




(B.7) 
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and therefore from Eq. (B.6) and (B.7) 


= Pkk + (Rfek AU)“^ - AU 


n 


<^kk - ^kk.xj“ + (R^^ AU)^^ 


(B.8) 


where ctnd the linearization Jacobians as defined in Eq. (B.6). 

Evaluation of these Jacobians will show that for constant transport coeffi- 
cients 

a 


kk 9x, kk 
k 


R,,_ = 0 


(B.9) 


and thus 


^^kk = air <®kk 


(B.IO) 


If the approximations outlined above are introduced into Eq. (B,3), we obtain 


AU 


n _ QAt 9 
1 + ^ I 9x 


+ ^ (Rjl AU)" 


6At 

1 +e 


-b" AU“+ |^(R22 AU) 


-C" AU“ + ^ (R33 AU) 
n-1 


^(AVj3+ AV,,) 


13' 


i 

i 




■*■ dz ‘^^ 32 * 


(Continued) 
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+ ^21+ ^22+ ^ 23 )” 


+ ^(-°+ ^31 + ^32+ ^33)”] 


+ T+1 


(B.U) 


Thus, for constant transport coefficients only the linearization Jacobians 

are needed in addition to the A, B and C Jacobians, 

Expanding and rearranging Eq. (B.ll), we obtain 


QAt fa n a a 


a „n 9 


+ ^c“ - -2-2 R 


3y 


2 22 


^ 1 

33 J 


* AU 


n 


= !f|fe <^^12 + ^^13)"’' + ^ <^^21 + ^^23 

^'^''31 ■*■ '^'^32* ] 


n-1 


+ ^(-^ + ^21 + ^22 + ^23)'' 


+ ^ (-G + V31 + V32+ ^33)"] 




(B.12) 


no 



Note the special notation used in writing the left hand side (LHS) of Eq. (B.12) 
which really must be considered an "operator," operating on AU^, and which 
is of the form 

EHS(12) = jl+a + b+c|* (B.13) 

This can be written in a spatially factored form 

LHS{12) = j{I + a) (I + b) (I + c)j * Au" 

= |(I + a + b + c) 

+ ab + ac + be + abcl * AU^ (B*14) 


if we note that ab, ac, be, and abc all are at least an order of magnitude (in At) 
smaller than a, b, c. Thus, without loss of accuracy. 


LHS(12) = 






1 + g ' 3x 


9x 


■ . 9At , 9 „n 9 

I + TT-r ( ntr S ; 

8x 


l + g ' 9y 




9At . 9 ^n 
+ g ^ dz ^ 


9z 


R 


R 


»>] 

«>]{ 


« AU 


n 


(B.15) 


Following Beam and Warming (Ref, B-1), in practice Eq, (B,15) is implemented 
by the sequence 


I + Ml /_9 .a“ - r" 
9x^ "y 




* AU = RHS(12) 


^ ^ 1 + ^ l9y ® 


ay 


R 


n\l 

/\ 


* AU = AU 


(B,16a) 


(B,16b) 


where RHS(12) means the right hand side of Eq. (B.12). 


Ill 



(B.16c) 


I + 


9At 

i+i 


(_L c" r“ \ 
az^ 33/ 


,n 


i= AU = AU 


* 


U 


n+1 


+ AU^ 


(B.16d) 


The remainder of the analysis follows that of Ref. B-1. 


112 



B.3 SUMMARY OF EQUATIONS FOR UNSTEADY ALGORITHM 

The vector of conserved variables, U, and the flux vectors of Eq. (B.l) 

are 




P 

pu 


m. 

pv 

— 

n 

pw 


q 

_PE_ 




Pu^ 


m 

pu + p 


m Vp + P 

puv 

= 

mn/p 

puw 


mq/p 

u(pE + p) 


(m/p) (r+p) 



pv 


n 


puv 


mn/p 

F = 

pv^+ p 

= 

n Vp + P 


pvw 


nq/p 


v(pE + p) 


(n/p) (r + p) 



pw 


~ q 


puw 


mq/p 

G = 

pvw 

=: 

nq/p 


pw + p 


q Vp + P 


w(pE + p) 


(q/p) (r + p) 
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where the pressure is given by the equation of state 


P 


P(Y-l) 


/ 2 ^ 2 ^ 2 
/ u -H V + w . . 


= ( 7 - 1 ) 


2 

m 


. 2 ^ 2 \ 

+ n + q \ 

2p j 


The viscous flux vectors are 


11 


0 

(2]Lt + X) u 


X 


pv 


X 


X 


/iw 

(2jLt + X) uu^ + + /Xww^ + kT^ 


V 


12 


0 

Xv 

> 

JLLU. 


y 


XUV + LLVU 

y y 


V 


13 


0 

Xw 

z 

0 

Xuw + /iwu 
z z 
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21 


Pv, 


X 


Xu 

X 

0 

puv. 


+ Xvu 


X 


J 
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In order to write the viscous flux vectors in terms of the conserved 
variables, the temperature gradients must be expressed in terms of the con- 
served variables. It is easily shown that 



Using this equation it can be shown that 


V 


11 


0 

(2fi+\) (pm^ - mp^) 

pp‘^ (pqj, - qp^) 

+ (1 - ^) [n(pn^ - np^) + q (pq^ - qp^)] 

+ ^ P(PJ^x ■ ''^xM 


V 


IZ 


0 

-2 

Xp (pn - np ) 

_2 y y 

MP (Pr^y - rnPy) 

0 

-3 -3 

AP m(pn^ - np^) + pp n(pm^ - mp^) 
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V 


13 


0 

Ap"^ (pq^, - 
0 

PP’^ (pm^ - mp^) 

Xp"^ m(pq^ - qp^) + pp'^ q(pm^ 


V 


21 



<P"x 

(pm^ 


np^) 

mp 


) 


0 

- 3 - 3 

pp m(pn - np ) + Ap n(pm 

X X. 


V 


22 


0 

PP’^ (pm^ - mp^) 

( 2 p+ X) p'^ (pn^ - npy) 

pp‘^ (pqy - qPy) 

PP'^ (2 + ^ “(P"y- "^Py) 

+ ( 1 - ^> [m(pm^ - mPy) + 

+ ^ P (PJ^y - ^Py)| 



0 

0 



(pq^ - '5Pz> 

(P"z - “Pz> 
n(pq^ - qP^) + PP 


q(pn 


-mpj 


- tnp ) 


q(pqy - 


- np^) 



0 


(P^x " 

0 

Ap"^ (Pr^x “ 

/xp"^ m(pq^ - qp^) + Ap ^q 


V 


32 


“O 

0 

pp’^ (pqy - qPy) 

^“^(PHy-nPy) 

jLip'-^ n{pqy - qPy) + Ap q(pn^ - np^) 


33 


0 

MP'^ (pm^ - mp^) 

PP’^ (P«2 ■ 

(2p+X) p'^ (pq^ - qp^) 

P + I 

+ (1 - i^) ['"(P'nz; - + "(P”z - 

+ ^ p(pr^ - rp^)j 
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The linearization Jacobians A, B and C are: 


A 


3E 

= '5u 


0 

v-1 , 2^ 2^ 2. 2 

(u + V + w ) - u 

-uv 

-uw 

u['}E-|(pE+p)] 


1 

* (y- 3)u 

V 

w 

i(PE+ p)-{y-i)u^ 


0 

- (y- i)v 

u 

0 

- (y -1) uv 


0 

- (y-l)w 

0 

u 

- (y-l)uw 


(y-1) 

0 

0 

yu 


® = -Su 




/ “ . 2 , 2 , 2 
(u + V + w ) ' V 


^[yE-|{pE+ p) 


0 

V 

(y-l)u 

0 

(y-l)uv 


- (y-3) V 
w 


i(PE+ p)-{y-l)v^ 


0 0 

0 0 

- (y-l)w (y-1) 

V 0 

- (y - 1 ) VW y V 


C 


9G 

= -5u 


0 

-uw 


-VW 

(u^ + 

w|yE-|(pE+ p)] 


0 

w 

0 

- {y-l)u 

- (y - 1 ) uw 


0 1 O' 

0 u 0 

w V 0 

- (y-l)v -(y-3)w (y-1) 

-{y-l)vw ^(PE + p) - (y-1) yw 
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The linearization Jacobians ^^22’ ^33 


R 


11 


av 


11 


au 


M 

p 


(2 + ^)u 


w 


(1 *^) (U^+ v^+ 


JL 

Pr 


E + (1+-) u 


-(2 + ^) 




0 

0 

-1 

0 




0 

0 

0 

-1 

(1 


0 

0 

0 

0 

Pr 


r 


R 


22 


av 


22 


9U. 


M 

P 


(2 + ^)v 


(1 -^) (U^+ v^+ w^) 


Pr 


E+ (l + ^)v‘ 


0 

-1 

0 

0 


-d~)u 


0 

0 

-(2 + ^) 

0 


0 

0 

0 

-1 

(1 -^) W 


0 

0 

0 

0 

.JL 

Pr 


R 


33 


av 


33 


au 


P 


u 


(l-fi 
JL 

Pr 


(2 + |)w 


y 

“) (u 


+ w^) 


+ ^ E+ (l + -^)w' 


0 

-1 

0 

0 


0 

0 

-1 

0 


0 

0 

0 

/•> A ^ y ^ 

(2 + ;i-rr)^ 


0 

0 

0 

0 

JL 

Pr 



b.4 three-dimensional steady parabolized algorithm 

A flow model which uses a spatial forward marching procedure in the 
principal direction of flow to obtain a solution cannot tolerate the upstream 
propagation of any flow phenomena. Such a model is obtained by deleting 
those viscous terms from the governing equations which contain gradients 
in the marching direction, and the resulting set of equations is termed 
"parabolized.” 

The three-dimensional, steady state compressible conservation equa- 
tions for parabolic flow (in the x-direction) are 

If + If + I? = ^ [v22(E- V ^ 

+ ^ b2<E.Ey) + 

where E is the vector of conserved variables, and F, G and V.^ are flux 
vectors , 


In complete analogy to the treatment of the unsteady problem, a general 
ized single-step spatial finite difference scheme for advancing the solution of 
Eq. (B.17) can be written as 


AE 


n 


QAx d . „n 

T+1 ^ ^ + 


Ax 9E^ 

1+^ ax 1+^ 


AE 


n-1 


where E^ = E(nAx, y, z) and AE^ = - E^ (terms of order Ax^ 

have been negelected for simplicity). 


and Ax 


(B.18) 

3 


Solving Eq.(B.17) for dE/dx and substituting into Eq.(B.18) yields 
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= 111 [ 4 ^ + ^^22 + ^^ 2 / 


+ ^ (-AG+ AV 32 + AV 3 / 


Ax 
1 + ^ 


^(-F+V22+ V23)” 




+ al (-=+^32+^33) 


n -1 


n 


(B.19) 


Again, AF , AG and AV-. are nonlinear functions of the conserved variables 
E. A linear equation with the same spatial accuracy as Eq. (B .19) can be ob- 
tained by expanding AF^, AG^ and AV^^j^ in a Taylor series while treating 
Avf*^ (i / j) explicitly. Accordingly, 


Af“ = (iff AE“ = A“ AE“ 

ag" = (lif ae" = b“ae“ 


(B.ZOa) 

(B.20b) 


where A and B are linearization Jacobians. 


While assuming that 


we can write 


AV?. = Avf.”^ 

ij ij ' ^ 


AV 


n 

kk 


'av 


n 


n 


kk 


9E 




• .-.n . _ n . -.-.n 
^kk^^ +^kk^^x. 


(B.21) 


(B.22) 
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since 


r" ak“ 

®kk3x, ^ 




we can rewrite 


^^kk = (Pkk - air «lk> + air <^lk 

k k 


(B.23) 


where and a-re the linearization Jacobians defined in Eq. fB.22). 

Using the viscous flux vectors V ^2 and can be shown again that, 

for constant transport coefficients. 


® r" = 0 


.n 


kk 9x, kk 
k 


which allows us to simplify Eq. (23) to 


(B.24) 


^^Ik = ^ («kk 


(B.25) 


Substituting Eqs, 


l^l-''*''22*''23l” 

■1 


6Ax 

1 + e 

Ax 

1 + 1 




+ T+i 


n-1 


(B.26) 
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From Eq, (B.26) it is concluded that in addition to the A and B Jacobians we 
only need (k=2,3) if constant transport coefficients are assumed. 


Expanding and rearranging Eq. (B.26) to combine all terms containing 


.n 


AE , we obtain 


I + 


0Ax 

1+e 


0Ax 

1+e 


2 2 

ay ^ ^2 ^22 ^ az ^ . 2 ^33 


ay 


dz 


* AE 


n 


_L Av“-1 + -9- Av“’^ 

9y 23 ^ dz '^32 


+ ^ [W ^ ^ ""22 ^23)“ + ^ (-C: + ^32 + ^33)“ 




n-1 


(B.27) 


In analogous fashion to the unsteady formulation and without loss of 
accuracy, we can write 


LHS(27) = 


T + /_L - -L- ^ 

' ^ 1 + 1 U ^ 9y2 «22, 


1 + 1 Uz g^2 33; 


* AE 


n 


(B.28) 


where LHS(27) means the left hand side of Eq. (B.27), which in practice is 
implemented by the sequence 


I + 


I + 


and, finally 


9Ax j 
1 + 1 1 

f—A^' 

[dy 

a^ 

8y^ 


* * 
AE 

0Ax 1 

(_L 


»»)' 

AE^ 

1 + e ' 


'az^ 


= 

-.-.n 
E + 

AE 


(B.29a) 


(B.29b) 


(B.29C) 


where RHS(27) means the right hand side of Eq. (B.27). 
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B .5 SUMMARY OF EQUATIONS FOR STEADY STATE ALGORITHM 

For Steady flow, the vector of conserved variables, E, and the flux 
vectors of Eq. (B. 17 ) are 


pu 

2 


^1 

pu + p 



puv 

= 

^3 

puw 



u(pE + p) 




pv 


E3/U 

puv 



2 

pv + p 

= 

^ 2 

^2 " EjU 

pvw 


E3 E^/ (E jU) 

v(pE + p) 


E3 E^/CE^u) 


pw 


p 

I 

puw 


^4 

pvw 

= 

E3 E^lEjU) 

2 

pw + p 


E2 + E^/(EjU) - EjU 

w(pE + p) 

U - 


E^ E5/(EjU) 


where u — f(E E^f E^, E^) as obtained by decoding the E vector for the 
primitive variables. It should be noted that although the flux vectors for the 
steady case are the same as for the unsteady case, their form in terms of the 
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conserved variables differs from that in the unsteady case. It is assumed 
that decoding the E vector was accomplished by obtaining u as a function of 
the E vector components, i.e,. 



and 

V = E3/Ej 
w = E^/E^ 

P = Ej/u 
p = - Ej u 

E = (E^ - up)/E^ 


It can be shown that in the decode procedure the (+) and the (-) sign apply to 
supersonic and to subsonic flow, respectively. As long as the flow is strictly 
supersonic or subsonic, choice of the sign should be no problem. It is ob- 
viously a problem for transonic flows, boundary layers and supersonic flows 
with imbedded subsonic pockets. 


In terms of the conserved variables, the heat conduction term becomes 


k 


9T 


JL A. / ^ 1 JL _L 

Pr dg \e J " Pr di 





Using this equation the viscous flux vectors V_ in Eq. (B.17) can be written 
as shown on the following page. 
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0 


22 


M 


S , ^ \ dK . 

E l au \ L 

9 y 

i=l 


-2 / ^^3 

iZfl + X ) ( Ej -^ 


- E 


OEj 

3 9y 


/e 


9E 


1 9y 
5 


4 . E 

^4 8y I 


, 5 , 3 . 8 E . 

<‘<■-*'"2 rar)TF 

i - 1 ' ^ 




i=l 

+ ^ + u ‘ ^ 


^ ^ ) E"^ E 


/ aE ^ 9 Ej^N 

3 (^1 “97 ' ^3 ~Wj 


+ -^ e" 

^ Pr 1 


-2 


+ M (1 - ^) Ej ^ ( Ej -^ - ^4 

9E, 


9E 


9E 


- E 


1 


9E 


- E, 


1 9y "^5 9y 


and 


23 


0 


-2 / 9^4 

(^1 9z ” ^4 9z 


-2 / ^^3 

me/(e 


9E 


1 az " ^3 9z 


/lE 


-3 / ^^3 

1 


9E 


E 


3 9z 


-3 / ^^4 

+ X Ej ^3 (^1 9z ’ ^4 9z , 


'J 
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wu ■uiHii mil 



and 
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The linearization Jacobians A and B are obtained by differentiating 
the flux vectors F and G with respect to the conserved variables. The 
result is 


> 

II 

> 

1 

> 

8u 

^9E. 

J 

13 


B = B1 . 

- Bl. 

9 u 

^9E. 

3 




where (9u/9Ej) is a column vector obtained by differentation. 



0 

0 

1 

0 

0 “ 


0 

0 

u 

0 

0 

Al. = - 
13 u 

-(u^+ v^) 

u 

2v 

0 

0 


-vw 

0 

w 

V 

0 


-v(E + ^) 

0 

(e + £) 

0 

V 



“ 1 

1 

1 

1 

1 


0 

0 

0 

0 

0 

P- 

2^ 2 
u + V 

2^ 2 
a + V 

U + V 

u + V 

2^ 2 
u + V 


V 

V 

V 

V 

V 


w 

w 

w 

w 

w 


(e + £) 

(E + £) 

(E + 

(E + £) 

(e + £) 
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a|Q. 



and 


BV. 





0 

0 

0 

1 

0” 


0 

0 

0 

u 

0 

1 


0 



0 

u 

“VW 

w 

V 


2 2 
-{u + w ) 

u 

0 

2w 

0 


-w(E + £•) 

0 

0 

(e + £) 

w 


1 

1 

1 

1 

1 

0 

0 

0 

0 

0 

V 

V 

V 

V 

V 

2^ 2 
u + w 

2 2 
u + w 

2^ 2 
u + w 

2^ 2 
u + w 

2^ 2 
u + w 

w 

w 

w 

w 

w 

(E + f) 

(E + f) 

(E + f) 

(E + f) 

(E + £) 


The components of the (8u/9E) column vector are given by 



y 

(y+i)E^ 



(±) ^ 


\Ai(h 



9u 

9^2 


(y+l)Ej (* ”Ej 



(±) 


TEj \^J 
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9u 


(±) 


(v-l)» 

yEj 



where 


3u _ i-r) (v-l)n 

9E5 - rEj 






,- 1/2 


The linearization Jacobians (h = 2, 3) are obtained by differentiating 

the viscous flux vectors with respect to the spatial derivatives of the 

conserved variables Ej^ (i= 1,5). The result is: 


R 


22 


= R 


22 




R 


33 


= R^3 + Q(l^) 


9E- 


■*^22 pu 


0 0 
0 0 


-(2 + -)v 0 


-w 0 


0 

0 

(2 + -) 
M 

0 


“I 


-d + ^)v" 0 

/I 1 \ / 2 , 2 

-^) (v + w ) 




' Pr' 


1 

Pr 
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Appendix C 

MULTILINEAR INTERPOLANTS FOR 
GIM MARCHING METHODS 


by 

John F. Stalnaker 



Appendix C 


C.l INTRODUCTION 

The following is a study to determine the nature of the multilinear 
weight functions which will generate certain implicit, spatial marching 
finite difference schemes within the GIM framework. The derivations are 
performed using rectangular two- and three-dimensional grids for sim- 
plicity of understanding. In these grids the local and Cartesian coordinates 
coincide; however, it should be realized that this is not always, perhaps 
seldom, the case and that the finite difference scheme generated by the 
GIM code is in terms of the local coordinates. 

With the above caveat aside let us proceed by setting in one place 
the notation to be used herein: 


a 

W 


(k) 
a(3 

f. g, h 


D 


E,F,G 


£ * 0 , 0 

a, b 


shape function for element point a 
weight function for element point p 

element difference operators for 77^ direction 

unique components of and 

respectively 

vectors of conserved variables 

Beam-Warming marching parameters 

finite difference parameters in the rj^ 113 
directions, respectively 

determinant of the Jacobian of transformation. 
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II fill n mini IT 


nmnmfiiri f i 


r 


Subscripts 


i,j,k assembled grid point indices for rj^* 

respectively. 

= 2-D 

- 1, ...» 27 3-D 

a, (3 element point indices = 1, . . . , 4 2-D; = 1, ... 8, 3-D 

In all the following rj j (and x in the case of the rectilinear grid) is 
assumed to be the marching direction. 


C.2 TWO-DIMENSIONAL BILINEAR WEIGHT FUNCTIONS 

For the two-dimensional case the shape functions are assumed to be 
the same bilinear shape functions now in the GIM code (Ref. C-1): 


S = ^ 
a 


(1 " 

(1 ’ r]^) 

^1 ^2 
(1 - Til) H2 


(C.l) 


The weight functions are assumed to have the form 


= ''po - Cp 


1 '^l ' + S3>ll 


(C.2) 


where the are to be determined. The nodal numbering system for the 
individual elements and the nine node box are shown in Fig. C-1. 


The two-dimensional element difference operators are 


D 


( 1 ) 

(3a 






8S„ „ 9S . 

^ 9y a dy 


dr] dn^ dr] 


(C.3) 
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1, J + 1 


1+1, j+ 



Fig. C-1 - Two-Dimensional Nodal Numbering System 

(a) Element System; (b) Nine-Node Rectangle 




(C.3) 

(Concl’d) 


(2) ^ 

pa 



W. 


as 

dr] 


a 


dx 

drjj^ 


as 


a 


dr] 


dx 

9t72 


Noting that for the rectilinear element 


ax 
ar? j 


= Ax, 


_§2. 

ar)2 


Ay, 


ax _ ay 


0 , 


Substituting Eqs. (C.l) and (C.2) into (C.3) results in the following element 
difference operators 


If = So - " Si ■ " ^2 S3] 

~ ' ■ °P4 ■ 12 Ax So ’ ^ Si ' ^ ^ ^ S 3 ] 

(C.4) 

^3 1 — 7)(^) _ _ 15^3) _ 1 r6c -2C -3C +C 1 

Ay “ S^ ■ Si ■ 12 Ay I® So ^ Si ^ ^32 ^ ^ 33 J 

^ = d(2) = - d( 2> = ^ [6 Cpo - 4 Cpj - 3 Cp2 H- 2 

The differential equation 

^ + ^ = 0 

3x 9y 

is modeled with a general spatial marching scheme (after Ref. C-2) 


+ ^ ^i+1, j+1 + <1 - ^i+1, j ^ ^i+1, j-l] 

+ ^1,3+ 1 + ^i,3-ll 


(C.5) 
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- [a F. , + (1 - 2a) F, , . + (a-1) F. , . 


(C.5) 

(Concl'd) 


= 0 


Examination of Eq. (C.5) reveals the significance of the difference parameters; 

a cross -plane {r} 2 ) difference parameter 

= 1 for forward cross -plane differences 

= 0 for backward cross -plane differences 

= ^ for centered cross -plane differences 

e marching (rjj) difference parameter 
= 0 for forward differences 

= -1 for backward differences 
= -i for centered differences 


0 weighting parameter for plane (i+T 


• th 


> 0 cross -plane differences are taken in the i 
and (i+1) st plane 

= 0 cross plane differences are not taken in the 

(i+1) st plane 


.th 


weighting parameter for plane (i-1) 

> 0 cross -plane differences are taken in the i' 
and (i-1) st plane 

= 0 cross -plane differences are not taken in the 

(i-1) st plane 


Note that for the present explicit differences in the elliptic GIM code 6 = ^ = 0. 


Assuming a form of the weight functions similar to that presently in 
the code, i.e., 

- ^ 


W 




(C.6a) 
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results in 

^(31 " s' ' ^(32^1 

%Z "" 3 " ^(32^/^^(31 “ ^132^5 

and 

^pi ■*■ ^|32 = Spl ■*■ ®p 2 

Assembling the elements as in Ref. C-1 in terms f ’ s and 

these to the difference coefficients in Eq. (C.5) and substituting 

• 

yields weight functions of the following form': 

(Pi ■ ’’P (^/^ ■ ^z'> 

iX 

^2 = ^(P2-”i)(2/3-'’2) 

J 

^3 = A" (P 2 ■ ” 1 ) (1/^ • 'l2> 

j 

W 4 = ^{Pi -rij) (1/3 -r,^) 

where 

= 36a(l+£ - 20 ) 

- 36 a (e - 2^) 


^Note that a' s 
similarity of 


and ( 3 's are substituted here for the d^^. in Eqs. (C 
these to the current GIM weight functions. 


(C.6b) 


j's, equating 
into Eqs. (C. 6 ) 


(C,7a) 


(C.7b) 


. 6 ) to show the 
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(C.7b) 

(Concl'd) 


= 36 (a-1) (e - 2<^) 


= 36 (a-1) (l + e-20) 


and 



(C.7c) 


These apply except in the case where 

e = 26 - 1; 6 0 

(e.g., Crank-Nicholson or ^ = 0, 6 = l/2). In these cases the weight functions 
do not maintain the same symmetry as the present weight functions. Values 
of the a's and p*s are available from the author. 

C.3 THREE-DIMENSIONAL TRILINEAR WEIGHT FUNCTIONS 

The procedure for the three-dimensional case is much the same as the 
two-dimensional problem. We assume the same trilinear shape functions now 
in the GIM code: 


(C.8) 
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a 


(l-r?j) (1-113) 

111 ^^"’^ 2 ^ 

T?! il2 
(l-i1l) ^2 
(1-rj j) ( 1 - 112 ) ^3 
r?l (1 -t?2) 113 

Hi H2 H3 
(1-Hl) 1I2 ^3 


I 


The weight functions are assumed to have the form: 

’ ^(3 "" “ ^(31 “ ^( 32 ’^2 " ^ jB 3 ^3 '*■ ^^2 

(C.9a) 

+ Cp5r]ir,3+ Cp^ r]2r73 - 77^172^3 

or 

^ ^^(31 “ ^1^ “ ^2^ ^^P3 ~ ^3^ (C.9b) 

The relations between the C’s and the d* s is obvious. The numbering system 
for the rectilinear element and the 27 node box are given in Fig. C-2. 


The full expressions for the element difference operators are given in 
Ref. C-3, For the rectilinear box, where 


At = Ax Ay Az, _ 

J drj 


ax 


= Ax, 


r^ rs ax. 

— ^ = Ay, = Az, - 0 i / j. 


1 


3t7 


Bin 


d-q. 


These operators become 

,( 1 ) 


pa 

= Ax 

pa 


Ay A^^y* drjj f dr,^ J" dr, ^ ^ 

O o o 

^ ^ as 

drii / ‘^‘’2 / 
o o o 

Ay ^ J* 


3 p ar]^ 


(C.IO) 


Substituting Eqs. (C.8) and (C.9a) into Eq. (C-10) results in four unique differ- 
ence operators for each direction 


fgj = Ax • D^y = -Ax . D<\>; ggj s Ay = -Ay • = Az 


d[V = -Az . 


P5 


pi 
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It can be seen from the above that 


" ®pi ®p2 ®p3 ®p4 ^pi **p2 ''p3 


and 


^P1 ^p2 “ ®pi ®p2’ ^p3 ■*■ ^p4 “ ®p3 ®p4 


^pl ■*■ ^p3 “ ^pl ■*■ ^p2' ^p2 ■*■ ^p4 ■ ^p3 ^p4 


®P1 ■'■ ®P3 = ^pl ^p4' ®^p2 ®p4 " ''pZ * ^p3 


Now, Eqs, (C.9b), (C.ll) and (C.12) can be combined to yield 


^B1 ^ fp3_. ^ ^ ^ ^ ^ 

^P2 ^p4 ’ ®p2 Sp4’ hp2 hp3 


-ipo = 




^Pl'''^p2'’'^p3'*'^p4 


* ^81'*' ^82* ' ^^63'*’ ^84* 


^P1''‘’^P2'*'*^P3'^’^P4 


^®pi '*' ®p3^ ■ 


*®pi ®p3^ ■ *®p2'*' ®p4* 


" 3 [^*®pi ■*■ ®p3^ ■ '®p2'*' ®p4* / 

'^P2 = i[^%l''‘'^p2> ■ %3'^^p4*P[%l'^^p2> ■ %3+^p4*] 

= T f2(f« , + ,) - (f„ , + 4^)]/f(f„ 1 + f« ,) - (fa , + f« .) 


P3 


pi p2' ' p3 p4'r I' pi p2' ' p3 p4‘ 


The differential equation 


^ M 4. ^ _ f. 

dx ^ dy dz " ^ 


is modeled by 


+ hp^ (C.lZa) 


(C.12b) 


(C.13a) 

P2'*' 


(C.13b) 
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s [(1 + - (1 + k + ^ E;;kJ 

+ ^ ^j«. k + (1 - ^l!k + (^ - 1) ^J!!. k] 

+ [- k + (1 - 2 a) f}^ k + (a - 1 ) fI,J 

- ^ [" i;}. k + (1 - 2a) + (a - 1) J (C.14) 

+ £ + (1 - 2b) 0^1 + (b - 1) 

+ [bo}^ + (1 - 2b) g}^ k + (b - 1) g;, k.J 


. A. 


i~i 


i-l 


,i-l 


^['^°j.k+l+ (l-2b)G.^j^+ (b-DG.^^.j] = 0 


where the difference parameters have the same significance as before and b 
is equivalent to a for differences in r] y 

Assembling the elements and equating coefficients does not lead to 
expressions for the weight functions in as straightforward a manner as in 
the two-dimensional problem. In order to resolve several ambiguities, the 
following considerations along with Eqs. (C,13a) were used: 

1. The weight functions should reduce to the form presently in the 
code for 0 = ^ = 0. 

2. The first four weight functions should readily reduce to the two- 
dimensional case. That is, the internal symmetries of the two- 
dimensional element difference operators should carry over to 
three dimensions. 

3. The weight functions derived here should be applicable to the 
elliptic solver with 0 = <^ = 0. Thus, the boundary terms must 
be differenced consistantly. 
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From this, the following weight functions result: 

cc 

Wj = ^ “17 j) (j “ 172) “ ^3) 

tJ 

ct 

= A~ (P 2 “ ^ 3* " 2^ ^ s’ “ 3^ 

J 

CL 

^ {{32 " 17 x) (s' “ 172) (s’ “ ^^3) 

J 

W4 = ^ ((3 X “ 17 J ) (^ - T7 2) (s’ ~ n 3) 

ij 

^5 = ^ (Pi -’ll) (§->’ 2 ) (I-I 3 ) 

^6 = ^(Pz'”!* (|''’2> 

J 

l) 

■^8 = ^(Pl''’l> (i'’^2> (i'l3> 

ij 


where 


a, = (a - b + b^) (1 + e - 20); a- 

i D 

= (a - b + b^) (e - 2^) ; 

= (a - 2b + b^) {e - 2^) ; 

= (a - 2b+ b^) {1 +e- 29); 


b(b- 1) (1 + e- 20) 
b(b- 1) (e- 20) 
(1-b)^ (e-20) 
(1-b)^ (l + e-20) 


and (3 X and are defined in Eq. {C.7c), 
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Appendix D 

VECTORIZED BLOCK TRIDIAGONAL EQUATION 
SOLVER FOR THE GIM/STAR CODE 


by 

S, J. Robertson 




Appendix D 


An attempt was made to develop a vectorized algorithm for solving 
large systems of linear equations of the form: 


L. 

X 


,U. , + M.U. + N. , 1 U. , , = D. 
1 1-1 1 1 1+1 1 + 1 1 


(D.l) 


where the L, M and N elements are 3x3, 4x4 or 5x5 matrix blocks, and 
the U and D elements are three-, four- or five -component column vectors. 
The subscript i in Eqo (D.l) corresponds to nodal points in a computational 
grid, and the three, four or five dimensionality of the matrix and vector 
elements depend on whether the system of equations are for a one-, two- or 
three-dimensional flow field problem (see Section 3), The system of linear 
equations represented by Eq. (D.l) forms a block tridiagonal system. 

A solution algorithm was sought that would make use of the parallel 
processing capability of the STAR- 100 vector computer. The Gauss -Seidel 
relaxation technique, based on an iterated solution of 


= (1 - wm:^ (L. 

1 ' ' i 11 




(D.2) 


where w is an over -relaxation factor, is the only technique which we could 
find that permits a straightforward use of vectorized computer programming. 
The inverse matrix in Eq. (D.2) is evaluated for all M. prior to entering 

the iteration loop. Since each matrix block is dimensioned only up to 5x5, 
the inverse can be evaluated by direct algebraic manipulation or by a Gauss 
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elimination technique. For the time being, we have coded only the algebraic 
inversion, since vector programming can be used in this method. 

Separate subroutines were programmed for 3x3, 4x4 and 5x5 block 
tridiagonal Gauss-Seidel equation solvers. These are listed in Tables D-1, 

D-2 and D-3 as subroutines EQSOLi3, EQSOL4 and EQSOL5, respectively. 

The argument list in these subroutines is (U, L, M,N, D, NODES, W, EPS, 
MAXI). The vector U is the solution vector which enters the subroutine as 
an initial or trial solution and returns as the updated or final solution. The 
matrices L, M and N and the vector D enter the subroutine as constants. 

The scalar NODES is the number of nodal points, W is the over -relaxation 
parameter, EPS is the error tolerance in the convergence test and MAXI is 
the maximum allowable iterations. The U and D vectors are doubly dimen- 
sioned, and the matrices L, M and N are triply dimensioned. The first sub- 
script of both vectors and matrices corresponds to the nodal point index. The 
second subscript of U and D corresponds to the vector components, and the 
second and third subscripts of L, M and N corresponds to the matrix elements. 

As of this writing, these subroutines have not been evaluated, except for 
some very simple test cases. They have not been applied to realistic fluid 
dynamics problems where their usefulness can be determined. 


154 


Table D-1 


LIST OF SUBROUTINE EQSOL3 


subroutine EUS0L3tU»L»lvi«N*D*NOD£S*t*fEP^iMAXI ) 
REAL L»N»^ 1 , N 

O I MENS I UN U ( 9 » 3 ) » L ( 9 , 3 • 3 > » M ( 9 * 3 » 3 ) ♦ M I ( 9 » 3 » 3 ) ' 
£ N(9t3*3)*D(9>3) >DET(9)*UP(90) »DIF(9»3) 
COmMON/XMI/mI 


MI 

( 

1 

« 1 

» 1 ibNUDE^ ) 

= 

,*l ( 1 

»2 


2i<NUDES ) * 

|V|( 1 

♦ 3 ♦ 3^NUUES ) 





3> 











ivl ( 1 ♦ 3 t Eit-NODES ) 

^ M ( 

1 

♦ 2 ♦ 3^ NUDES ) 

MI 

< 

1 

O 

f liNODEU ) 

= 

-M ( 

1 

1 

2 

, lipNODES) * 

iM ( 

1 O » 3*NUDES ) 
















+1*1 ( 

1 ♦S’ lib.NUDE'Sl'X- 

|V| ( 

1 


2 ♦ 3^NUUES ) 

|V|I 

( 

1 

* 3 

» 1 iNUDEU } 

= 

Mt 1 

« 

2 


1 iNUDE3 ) * 

M( 1 

♦ 3 ♦ 2^I''JUDE^ ) 
















I*. ( 

1 ’3» I^InUDEo)* 

,vi ( 

1 


2 ♦ 2i^NUDcLD ) 

MI 

( 

1 

« 1 

» 2 ^ N U D E 3 ) 

= 

-.Vi ( 

1 

« 

1 

, 23NODE3 ) * 

1*. ( 

1 O ♦ S^NUDEo ) 





S 











+M < 

1 ♦ 3 1 EiNuDES ) -X- 

M ( 

1 


1 » s^nudes ) 

MI 

{ 

1 

*2 

♦ 2ibN0DE3 ) 

= 

M ( 1 

n 

1 

« 

1 SNODES ) * 

M( 1 

♦3’3aNOD£S ) 





* 











M ( 

1 ♦ 3 * 1 •i’NuOES ) ^ 

M< 

1 


1 ♦ 3^N0DES ) 

T''! 1 

( 

1 

< 3 

♦2^NOdE3 ) 

= 

-.Vl ( 1 

« 

1 

« 1 i.NUDE3 ) * 

|V1 { 

1 ♦ 3 ♦ 2i*NUDE3 ) 





S 











+ |V| ( 

1 O ♦ 1 snode^ ) 

M ( 

1 


1 ♦ 2^NUDES ) 

|V| 1 

( 

1 

♦ 1 

♦ 3^NUD E 3 ) 

= 

,*1 ( 1 

n 

1 

% 

2 io NUDES ) * 

|V. ( 1 

♦ 2 * 3^NUDE3 ) 





S 











— (V'l ( 

1 « 2 ♦ 2itNuDEc2. ) -X- 

j*. ( 

1 


1 *3^ IN DDES) 

M I 

( 

1 

O 

» 3^imUdE3 ) 

= 

-(V. { 

1 


1 

« 1 ibNUDE3 ) * 

fv, ( 

1 ♦ 2 ♦ 3+NuDES ) 
















+ .V. ( 

1 ♦ 2 ♦ 1 ^NuDEo 1 

| V '| ( 

1 


1 ♦ 3^NuDEu ) 

MI 

( 

1 

« 3 

* 3^NOdE3 ) 

= 

M ( 1 


1 

9 

1 ANODES ) -X- 

|V-> ( 1 

♦ 2 ♦ 2it>NUDES ) 





ib 











.VI ( 

1 ♦ 2 ♦ 1 ^t'v'^DES ) "X" 

.V|( 

1 


1 ’E^iNvUEU) 


DET ( 1 SNUUES ) =1*1 ( 1 » 1 ^ 1 ibisjUDEU ) I ( 1 » 1 » 1 ilMuDES ) 

i +)^l ( 1 , 1 » 2ii\ODES ) 1 ( 1 ♦ 2» libNUUEis ) 

S +i>l( 1 , 1 » 3i>lMODES ) I (1.3* ISNOUES) 

DO 50 1 = 1 » 3 

DO 50 J=l»3 

M I ( 1 » I , Ja.NUDEU ) = ,^U ( 1 . 1 « USN0DE3 ) /DET ( la^NUDE:^ ) 

EO continue 

NMi =N0DES- 1 
NM2=N0DES-2 
ITeR=0 
10 continue 
I TeR= I TER+ 1 
DO 100 1=10 
UP( 1 * I ) = ( 1 .-W ) -x-u ( 1 O ) 

UP ( 1 » I ) =UP { 1 » I ) -vV*N I ( 1 » I « 1 ) * (N ( 1 * 1 * 1 ) *U ( 2 ’ 1 ) +N C 1 ♦ 1 *2 ) *U ( 2 * 2 ) 
+N ( 1 ♦ 1 '3 ) *U ( 2 * 3 ) ) 
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Table D-1 (Concluded) 

UP ( 1 » I ) =UP ( 1 » I ) *U (2»1)+N<1’2’2 ) *U i 2 ♦ 2 > 

S +N( 1 i2*3) »-U( 2*3) ) 

UP( 1*1) =UP ( 1 * I ) “ /J-Jt-M I(i»I*3)*(iM( 1*3*1) *u (2*1) +N ( 1*3*2 ) *U (2*2) 

S +N ( 1 * 3 * 3 ) *U ( 2 * 3 ) ) 

UP( 1 * I ) =UP( 1*1) (iv|i ( 1 * I * 1 ) •?«■□( 1 * 1 ) +|Vii ( 1 , 1 ,2) *D( 1 *2) + 
ip M I ( 1 * I * 3 ) -«-D ( 1 * 3 ) ) 

UP { 2 * I SN;<^2 ) = ( 1 • - J ) *U ( 2 * I S>iNiM2 ) 

UP (2*1 ibNM2 ) =uP ( 2 * I ipiNM2 ) I ( 2 * I * 1 ^N>i2 ) '^' ( L ( 2 * 1 * 1 -pNi''i 2 ) -^^-U ( 1 ♦ 1 4>lMl''i3 ) 

ifc +L ( 2 * 1 * 2ipiNM2 ) ( 1 * 2ifai'JM2 ) +L { 2* 1 * 3ii>Nf''i2 ) *U ( 1 * 3i^Ni^i2 ) > 

UP ( 2 » I iNM2 ) =UP ( 2 * I ipNiVi2 ) -»A/*c-i I ( 2 * I * 1 i»i\M2 ) *- ( N ( 2 * 1 « 1 ^iMM2 ) *U ( 3 * 1 a.t^M2 > 
i +N ( 2 * 1 * 2ii*N|Vi2 ) * J ( 3 * 2iHMM2 > +N( 2 * 1 * 3ipNl^i2 ) *U ( 3 * 3ii>Nlvi2 ) ) 

UP (2*1 iNi'^2 ) =0^ (2*1 ii>N|vi2 ) - W-i^-|V| 1(2*1* 2ipNKi2 )->'«■( L ( 2 * 2 * 1 ipNl'-'l2 ) *U ( 1 * 1 iNM2 ) 
ifc +L ( 2*2*2ipNM2 ) ( 1 * 2S>NM2 )+L( 2*2 * 3ii^NKi2 ) *U ( 1 * 3^Nlvi2) ) 

UP( 2 » I SNrl2 ) =UH (2*1 ipNM2 ) -W/*iv; 1(2*1* 2il>NiM2 ) * ( N ( 2 * 2 * 1 i>NjVi2 ) *U ( 3 * 1 iNM2 ) 
ib +N ( 2 * 2 *2ii>NM2 ) J ( 3 * 2ipNM2 ) +N ( 2 * 2 * 3ipNivi2 ) -«-U ( 3 * 3^-Ni''l2 ) ) 

UP (2*1 iiMM2 ) =UP ( 2 » I iiNtjvi2 ) 1(2*1* 3ii>NPl2 ) * ( i_ ( 2 * 3 * 1 iplslM2 ) -^^U (1*1 itN|vi2 ) 

^ +L ( 2 * 3 * 2i»»NM2 ) *0 ( 1 * 2^bi\M2 ) +L ( 2 * 3 * 3ipNlVi2 ) *U ( 1 * 3ifcf'itvi2 ) ) 

UP (2*1 S>iNi'/l2 ) =UP { 2 » 1 ^NM2 ) *“ 1(2*1* 3ipiNi''i2 ) * ( i'n ( 2 ’ 3 * 1 -p|^ii''l2 ) "^'U (3*1 ibN^i2 ) 

* +N ( 2 * 3 * 2i*>Nfvi2 ) ( 3 * 2iiNM2 ) +N ( 2 * 3 * 3ipNJ'^.2 ) *U ( 3 * 3ifcNM2 ) ) 

UP ( 2 * I ipNiVl2 ) =UP ( 2 * I i.Nfvi2 ) +iAi* ( M 1 ( 2 * I * 1 SNfVi2 ) -«-D ( 2 * 1 iNi'^iZ ) +M 1(2*1* 2it'NiVt2 ) 
^ *D ( 2 * 2i»NM2 ) +1''! I i 2 ’ I * 3i'>NR2 ) *D ( 2 » 3it'i'^M2 ) ) 

UP (NODES* 1 ) = { l.~',V)*U( NUDES* I ) 

UP( NODES* I ) =UP ( imODES* I )-'w*NI (NODES* I * l ) (L( NODES* 1 * 1 ) *U ( NN I * 1 ) 
i +L<N0DES* 1 *2) *'j(NIVii *2)+L(N0DES* 1 *3)-«-U(NNl *3) ) 

UP( NODES* I ) =UP ( (MODES* 1 ) -»Ar)(-,viI (NODES * I * 2) * (L ( NODES* 2* 1 )-«-U (NMl * l ) 

i +L ( NODES * 2 * 2 ) *0 ( NIV| 1 * 2 ) +L ( NODES * 2 * 3 ) *U ( NtM 1*3)) 

UP ( NODES* I ) =UP ( imODES* I ) - 1 ( N0DE^« I * 3 ) ‘»'^ ( L ( NUDES* 3* 1 ) *U ( NN 1*1) 

S +L ( NODE-^ * 3 * 2 ) ^ Ni'^i 1*2) +L ( NODES *3*3) *0 ( InIM 1*3)) 

UP ( NODES* I ) =UP ( .modes* I ) +** (Ml (NODES* I * 1 ) *D ( NUDES* 1 ) +M I ( NODES* 1*2) 

S -«-D ( NODES * 2 ) +i^l I ( (nIODES * I * 3 ) *D ( NODEo * 3 ) ) 

100 continue 

DO 350 1=1*3 

D If ( 1*1 a>NC OeS ) = OP ( 1 * ISNODES ) -U ( 1 * I SNODES ) 

350 continue 
PMS=0 • 

DEL=0. 

DO 360 1=1*3 

DO 360 0=1 *nODEo 

Di:lL = DEL+D 1F(J*I )'i^DIF(0* I ) 

RMS = PMS+UP ( J * I ) ^^UP ( J * I ) 

3o0 continue 

DEL=SOPT ( del ) 

I-<|V|S = S0PT ( Pi''i3 ) 

TEST=DEL/PMs 
DO 400 1=1*3 

U ( 1 * I SNODES ) =0P ( 1 * I il^NODES ) 

400 continue 

1F( TEST. LtL-EPS. OR. I TEP.GE.NAXI ) RETURN 
Gu TO 10 

end 
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Table D-2 

LISTING OF SUBROUTINE EQSOL4 


ijUB^^CJUT 1 Nh. L,^ ( o » L » M 1 N » U » NOUt. » EPS » I'^iAX I ) 

P£AL 

D 1 MENS I UN U(9»4) >L(y»4*4) »|vl ( 9 * 4 * 4 ) *ivii ( 9 * 4 * 4 ) * 
i iM(9«4'4)*D(9>4) »DET (9)«UP(9»4) <DIF(9»4) 

CUMMON/XM I /M I 

M I ( 1 -• 1 * 1 j3 1>1Ud£^ ) =i'1 ( 1 * 2 »2^NUD£i3 ) * ( I'l < 1 » 3 * 3^N0DE^ ) ( 1 » 4 »4^fcNODES ) 

S - |V| ( 1 * 4 * ai.NODE'^ ) ■^1’'' < 1 ’ 3 ♦ 4 ^NUD£S ) ) -M ( 1 * 3 » 2 ^ 1 '^OUES ) * ( ( 1 • 2 ’ 3it.\ODEi> ) 

S *|V 1 ( 1 f 4 * 4i>NOD ti-u ) ( 1 ♦ 4 * 3^NODE'^ ) *i'’l (1*2* 4 SNuUEo ) ) +M ( 1 * 4 * 2 i^NODES ) 

S * (!''!( 1 ♦2«3ibN'JUE-^)'*^N( 1 * 3 ♦ 455 NODES ) -M ( 1 * 3 * 3iNSDE::i ) ( 1 » 2 » 4 !t>NODES ) ) 

( 1 *2* 1^NvJDE!'^)=“’I''i^ 1 *2’ 1 ^NuDE^ ) ^ 1 *3’3-*Ji'^'^DE— 1 »4*4 ibi"^UDE ) 

i -M ( 1 * 4 »3^NQDES ) ^«-Ki < 1 ♦ 3 * 4 it>iMUuES ) ) +iv, ( 1 ^ 3 » 1 SISiuDL^ ) * ( Mi ( 1*2* 3 *NUDES ) 

i *;Vl ( 1 * 4 » 4i>f'JODES ) -I'-i ( 1 » 4 » 3iNODES ) *MI ( 1 » 2 » 4 ^NODE^ ) ) -i-i ( 1 » 4 » 1 iNODLo ) 

it * ( Ml ( 1 » 2 ♦ S^ijNuuE-J ) ( 1 * 3 » 4iuNuDE'iJ ) ~mi ( 1 » 3 * S-dNoDE-^ ) *i’'i ( 1 ' 2 ’ 4i>NouE'^ ) ) 

i''! I ( 1 » 3 » 1 SNOdE^p ) =Ml ( 1 » 2 ' 1 ibNODES ) * ( i^i ( 1 * 3 » 2itNU^Ec> ) *Ml ( 1*4 ♦ 4is>NuDES ) 

S "■ivi( 1 *4»2 ^nQL)ES) ^mi( 1 *3’4iblMs_>UES) ) — Mi( 1 ^3* 1 *2«2ibNv^3tz.E) 

ib *jVi ( 1 * 4 1 44 »NQDt.C 3 ) -Ml ( 1 « 4 ♦ 2^NUDE^ ) *Mi ( 1 » 2 » 4 ^bNUDEE ) ) +iv| (. 1 * 4*1 SNuDE^ ) 

It ( Ml (1*2’ 2i»NULjE-> ) *1^1 ( 1 ’ 3 ’ 4^NODES ) -mi (1*3* 2itlNuDE-i> ) *mI ( 1 ’ 2 ’ 4 S 1 MUDES ) ) 

Ml I ( 1 ’ 4 ’ 1 ibNUDE^ ) = -Mi (1*2’ 1 icNuOES ) * ( i'^l C 1 * 3 * 2ibNoDE^ ) ( 1 * 4 * 3 itNUDEi 3 ) 

i>-Ml ( 1 * 4 * a-bNUDE-^ ) ^t«'i ( 1 ’ 3 ’ 3^NUD£i::' ) ) +i"i ( 1 ’ 3 ’ 1 ^b^NwUE^ ) * ( ,^i ( 1 ♦ 2 ’ 2-bi'J'-'QE^ ) 

it -K-|V| ( 1 * 4 * 3ibNGDE^ ) -Ml ( 1 ’ 4 ’ 2^NUQEP ) ^Mi < 1 ’ 2 ’ 3ibNuQE^ ) ) “Ml ( 1 * 4*1 itNuDEib ) 

S *(Mi( 1 * 2 ’ 2 ^NJUE'P)‘''^Ml( 1 *3*3ibNODE'ii')— i''l( 1 ’3’2ibNUDE^)*Mi( 1 *2’3^NJDE^> ) 
Mil ( 1 * 1 *2ibNUQEo) = — Ml( 1 * 1 * 2-bNODES ) ( I'-i ( 1 ♦3’3ibN'^DES)'ii'i''i( 1 *4*4^I^1->DEe) 

i' — M ( 1 ’ 4 ’ S^i'^OD EO ) -«-Mi ( 1 * 3 * 4ib|\JUiJEii^ ) )+!''! ( 1 ♦3’2-bNUUc.'b)-i<-(i''i( i * l * 3ibiNiODE-> ) 

it) *m( 1 *4*4i^>N0DEo) — i^i( 1 *4’ 3^NOUEiP ) *Mi ( 1 * 1 * 4ibt\iuut.o ) )— i''i( l *4*2iti'^PEE^) 

ib * ( i»i ( 1 * 1 * 3ibN JGE-J ) *mi ( 1 » 3 * 4-t NODE‘S ) -l'’i ( 1 » 3 ’ 3itiNUDEE' ) ^Mi ( 1 ’ 1 ’ 4 i>iN^^DEE) ) ) 

Mil ( 1 « 2 'E*^ -^DEo) = ,'|( 1 * 1 * 1 ojinUDES ) * ( Mi { 1 ♦ 3 ♦ 3 ) *Mi ( 1 * 4 * 4 ^NuijEe) 

it — M ( 1 * 4 * 3^NQl)eo ) ^Mi ( 1 * 3 ’ 4^NODE^ ) ) —Ml ( 1 * 3 ’ 1 -bNtjEE-b ) -«- ( Mi ( 1 * 1 * 3:»iNE0Eib ) 

it *m 1 ( 1 * 4 ’ 4ibN0D E ) — i"^! ( 1 » 4 ’ 3 it Nudes )*i»i(1»1* 4 ibNUDE-^ ) )+M|( 1 *4’ litNUDEU) 

S * ( Ml ( 1 * 1 * 3ibN'JDE'^ ) *Mi ( 1 ’ 3 ’ 4ibNUDES ) -Ml ( 1 » 3 ’ 3 ibi\UDEE ) -5t-|V| ( 1 * 1 * 4ibNuDEu ) ) 

Ml 1 ( I * 3 * 2 ^NudEu ) = — I'l ( 1 * 1*1 ibNuDEu )*(mi( 1 * 3 * 2 ^NUDEu ) -«-mi ( 1 * 4 * 4 iiiNuDEU ) 
it -M ( 1 * 4 ’2ib'NiODES ) -«-Mi (1*3* 4ibNODES ) ) +M| (1*3*1 iNUDEu ) * ( M ( 1 * 1 * 2^N0DEi=> ) 
it *M ( 1 * 4 * 4ibiNOD Es 1 —Ml ( 1 * 4 * 2ibNuDES )*mH 1*1*4 itN(JDEu ) ) —mi ( 1 * 4*1 sNUDEej ) 
it -K- ( Ml ( 1 * 1 * 2itNJuEE ) ( 1 *3 ’ 4it NODES ) -Mi ( 1 * 3 * 2 iti\UDEu ) *,v| ( 1*1 * 4itN0DES ) ) 

Ml 1 ( 1*4’ 2ibNDoEu ) =ri ( 1 * 1 ♦ 1 itNODEu ) * ( Ml ( 1 * 3 * 2 ibNUDEU ) ( 1 * 4 * 3itNODEu ) 

i -Ml (1*4 * 2ibiM0DES ) *Mi (1*3* 3ibNUDES ) ) -Mi (1*3*1 ibNUDES ) * ( MU 1 * 1 * 2iti\ODES ) 

43 *M1 ( 1*4 *3itiNjoDEu ) -Ml (1*4* 2ibNuDEU ) ^-i''i (1*1* Sit-NUDEE* ) ) +mi ( 1 * 4*1 ^jNODeu ) 
i) * ( Ml ( 1 • 1 * 2itNODE'iP ) *bi ( 1 * 3 * 3i»N(MDES ) -Ml (1*3* 2 ibNUDEE ) ^eMi ( 1 * 1 * 3 itNODES ) ) 

Mi I ( 1 * 1 * 3^NUdEu ) =m '1 ( 1 * 1 * 2 ibNuDEE> ) ■»• ( Ml ( 1 * 2 * 3ibNuDEi:j ) *Mi ( 1 * 4 * 4iNODEU ) 

it — Ml ( 1 * 4 ♦ JibiNOD eS ) ^Ml ( 1 * 2 * 4 ibiNiuDEU ) )— mi( 1 *2*2^i'^ODE'J)'^(MI( 1 * 1 *34jN>wiDEe) 

i *M1 (1*4* 4itNQDES ) -Ml ( 1 * 4 * 3 ^NODES ) *Mi (1*1 * 4itNODES ) ) +mi ( 1 * 4 * 2 itNuDEu ) 
it * ( |M 1 ( 1 ♦ 1 * 3itN0oE'P ) *Mi ( 1 * 2 * 4ibNODEu ) —Mi (1*2* 3ibNuDEu ) -M-ivj ( 1 * 1 * 4itN0DES ) ) 

M I ( 1 * 2 * 3 -tN 0 DE^ ) =-M ( 1 * 1 * litNODE^ ) * ( Ml ( 1 * 2 * 3ibN0DEU ) *|v| (1*4* 4SNODES ) 
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Table D-2 (Continued) 


^ -M ^ 1 «4 ’ 3-i-NODLS ) ( 1 ♦ 2» 4*NUDEij ) ) +M ^ 1 * 2 » 1 ifcNODE^ ) * ( M ( j » i , 3:i>|Nj(jQES ) 

i. 1 ♦4»4^N0iji£S ) -I'H 1 » 4 » 3itN0DE3 ) -K-|V| ( j , j »4ibi\0lJEo) ) -Ki ( 1 »4» ISNUDES) 

4a * ( I'i ( 1 ♦ 1 » 3^N3iJL'i5 ) *i''i (1*2’ 4^N0DE3 ) “ri ( 1 » 2 » S^I'^jODE^ ) ( 1 * 1 « 4ibN0DES ) ) 

MI ( 1 *3»3itNUDE^ ) =i'i( 1 * 1 * lipNUOE^ ) (M( 1 • 2 * 2^NOL)E>:3 ) *M ( 1 * 4*4iNODE^ ) 

— M (1*4* 2^NO0 Lo ) ( 1 * 2 * 4i-iM0DEo ) ) — |Vi (1*2*1 4al\UDL— ■ ) * ( I'l { i ♦ i t 2^M0 l)ES ) 

i *|vi ( 1 » 4 * 4iM0UL3 ) (1*4* 2iN0UE-j ) (1*1* 4^iNuuE^ ) ) +l’>'i (1*4*1 iNQDE3 ) 

is * C i''l ( 1 ♦ 1 *.2 4aN0L^E-::3 ) '^M (1*2* 4i>NODES ) -iM (1*2* 2i'NUDL3 ) ^i>i (1*1* 4iDNuULS ) ) 
i'i 1(1*4* 3itN^DE-i ) =-)'* (1*1* lii>N0DE3 ) * ( K ( 1 * 2 * 2isNODE3 ) *|V, ( 1 , 4 , 3^i\lODEi3 ) 

^ ~M( 1 *4*2^NQ[)E3) ^i'^i( 1 *2*3it*M00E3) )+i''i( 1 *2* liiaI\iuOEiJ)'^(l'''i( 1 * 1 *2it’NUDEij) 
ill ( i * 4 * 3^M0L)Lo ) — ( 1 * 4 * 2i*^NUDEi^ ) ( 1 * 1 * 3^M(jDEiiJ ) ) ~i''t (1*4*1 saNUDEi) ) 

iL * ( I'^i C 1 » 1 * 2 iJ E 3 ) *r'^i ( 1 * 2 * 3^ NODE 3 )— i''i( 1 *2*20f\UDL— )*i''i( 1 * 1 *3it>NODE'is) ) 

M I ( 1 * 1 * 4i>NODE3 ) =~M (1*1* 2ifNODE3 ) -K- ( j-i ( 1 * 2 * 33Nv^DE3 ) *,'/i (1*3* 4SNODEis ) 

E -.VI ( 1 * 3 » 3ON0DE3 ) -!(-Ni ( 1 * 2 * 4^N0DE3 ) ) +iv| (1*2* 23N0DE^ ) ■«■ ( IW ( 1 * 1 * 34 -NUDE 3 ) 

it *|Vi ( 1 ♦ 3 *40NODES ) -M ( 1 * 3 * 33N0DE3 ) -K-M ( 1*1 » 4isNUDE3 ) ) -ivi ( 1 * 3 , 2 *NUDE 3 ) 

in * ( i'i (1*1* 3it>N0 0£^ ) *)'“! ( 1 * 2 * 4^NODEis ) -|Vj ( 1 * 2 * 3^NuDEis ) -x-.v, ( 1 * 1 » 4itNODE3 ) ) 
M 1 ( 1 * 2 * 44aNUD£o ) =i‘i ( 1 * 1 * 1 itNODE'-> ) ■^'" ( 1^1 ( 1 * 2 * 3isNOUE3 ) ( 1*3 *4isNUDE;:j ) 

3 -M ( 1 * 3 * O^l'^ODE^ ) -«-|V) ( 1 * 2 * 4 ii'NOuEis ) ) -|V 1 ( 1*2*1 3NUDti-i^ )•«<-( M ( 1 * 1 * 34*iNODEi^ ) 

it' "W’|V|( 1 *3*4itNODL.S) ~M( 1 »3*3isiNUUE3)"'<'i''^l( 1 * l *4i*^NUDE3) )+i'^l( 1 *3* litNODEO) 

ib -'i- ( i*'i ( 1 * 1 * 3^NOlJL'-> ) ■" *’> (1*2* 4^NODEis ) — i"'! (1*2* 3^NODL^ ) (1*1* 4d>iMuUEi3 ) ) 

I ( 1 * 3 * 4iNODE3 ) — — M (1*1*1 itNUUEE ) ( M ( 1 * 2 * 23fNODE3 ) *i-'i (1*3* 4itiNjvJDLi> ) 

3 -M ( 1 * 3 * E^iMQUcS ) ^I'-i (1*2* 4isNxjDE3 } ) +iv, ( 1*2*1 ^NODE-i^ ) * ( ( 1*1* 2i^>N0DEis ) 

ip ( 1 * 3 *431301)03 ) -M (1*3* 23N0DE3 ) (1*1 * 43^JODE3 ) ) -)M ( 1 * 3 * 1 itNGDEo ) 

3 ■)(■ ( M ( 1 * 1 * 23N0L;£o ) (1*2* 43NUDE3 ) -|vi (1*2* 23KODEo ) -)i-|vi (1*1 » 43alMUDE3 ) ) 

13 I ( 1 * 4 * 40iN0DES ) =M (1*1*1 SNODES ) * ( M ( 1 * 2 * 24.N(JDE3 ) *M (1*3* 3itNODES ) 

it -M ( 1 * 3 * 2 ^NQDc:S ) -)(-M (1*2* 3-NUDE3 ) ) -ivi ( 1*2*1 itNuDLis ) -)i ( (3 ( 1 * 1 * 2 -pN0DE3 ) 

it *;M (1*3* 3itNCDE3 ) -13 (1*3* 2t>NUD£3 ) (1*1* 3itNUDE^ ) ) +,3 (1*3*1 isNODEtj ) 

* ( i»i ( 1 * 1 * 2 i^N^i^£o ) ( 1 * 2 * 3^1'jL'DEib ) “I'I ( 1 * 2 * 2-i^i^3iDE0) ) *i''i ( 1 * 1 * oaaiNiuDEiis ) ) 

DET ( 1 iNOOEE )=i-i{ i* 1 * 1 t>iN.0'DEE ) I ( 1 * 1 * 1 ifa(\UDE3 ) 

it +i3 ( 1 > 1 * 2 ^InODEE ) *i'vi I ( 1 ♦ 2 * 1 ONODEE ) 

i { 1 » 1 * 3iti\0DEE ) 1 ( 1 * 3 * 1 itNODEE ) 

it +M ( 1 > 1 * 4ii’MODES ) -)«-l'3 I ( 1 * 4 * 1 itNODEE ) 

Ni''i 1 =NODEE— 1 
Ni32 = N0UEE-2 
DO 50 1=1*4 

DO 50 0=1*4 

MI ( 1 ♦ I *04>NOq£E)=,'|I ( 1 * I * Jit NODES ) /DET ( 1 SNOuES ) 

50 COnTINOE 

iten=o 

10 Cunt I NOE 

I TEN= I TEN+ 1 

DO 100 1=1*4 

OP ( 1 * 1 ) = ( 1 . - tV ) -«-U ( 1 * I ) 

OP( 1 * I ) =OP (1*1 ) - I ( 1 * 1*1 ) * ( N ( 2 * 1 » 1 ) *0 ( 2 * 1 ) +N ( 2* 1 * 2 )*U ( 2 * 2 ) 

it +N ( 2 * 1 * 3 ) -)(U ( 2 > 3 ) +1'^ ( 2 * 1 * 4 ) *U ( 2 * 4 ) ) 

OP ( 1 * I ) =OP (1*1) - vJ-Jt-M I ( 1 * I *2 ) * ( IN ( 2 * 2 * 1 ) *0 ( 2 * 1 ) +N ( 2 * 2* 2 ) *U ( 2 * 2 ) 
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Table D-2 (Continued) 

S +N(2»2»3)*U(2>3)+N(2^2*4)*U(2*4)) 

UP( I * I ) =UP < 1 * 1 ) II I*3)*(tNl2’3*l)*U{2’l)+N(2»3»2)-^U(2»2) 

S +N(2»3»3)*U{2*3)+N(2»3*4)*U(2*4)) 

<JP ( 1 < I ) =PP t 1 * I ) I ( 1 * I * 4 ) * ( IM( 2 » 4 * 1 ) *U ( 2 » 1 ) +N ( 2* 4 » 2 ) -^ U ( 2 ’ 2 ) 

* +N(2»4*3)^U(2»3)+N (2*4^4) *U(2»4) ) 

UP ( 1 5 I ) =UP ( 1 > I ) ( M I ( 1 » I » 1 ) ( 1 * 1 ) +M I ( 1 » I * 2 ) ( 1 » 2 ) + 

S. M I ( 1 t I » 3 ) *D ( 1 » 3 ) I(1*I»4)*D(1*4)) 

UP C 2 « I iiiMM2 ) = ( 1 . - .fn *U ( 2 » I ii>NM2 ) 

UP ( 2 < I ) =UP ( 2 » 1 i±>Nr'l2 ) -W-«-tvj I ( 2 * I » 1 ^Ni''i2 ) * < l. ( 1 * 1 < 1 ipNi''i 2 ) *U ( 1 * 1 a»iNH 2 ) 

^ +L < 1 » 1 ’ 2^N)V.2 ) '^U I 1 * 2^I'4M2 ) +U ( 1 » 1 » 3-t'N.>'i2 ) *U ( 1 » 3a>jNKi2 ) +i_ ( 1 » 1 » 4i^l‘'l2 ) * 

i£ U ( 1 ♦ 4SNl''l2 ) ) 

'-«P { 2 » I ii'^)^l 2 ) =uP ( 2 » I ^NM 2 ) - wS-M I ( 2 » I < 1 ^inM 2 ) * ^ N ( 3 i l » i 4 jNivi 2 ) *U ( 3 * 1 ^Nr>'i 2 } 

^ +N < 3 » 1 *2^Nm 2 ) ( 3 ’ 2ii^NN2 ) +N( 3 * 1 » 3ibNM2 ) *U ( 3 » 3S"NM2 ) +N ( 3 » 1 *4^NiVi2 ) * 

S U ( 3 ♦ 4it>NM2 ) ) 

UP ( 2 ♦ I St>N|Vi 2 ) =UP { 2 ♦ I ifaNiVi 2 ) -W-X-[vi I ( 2 ’ I ’ 2ii>NPi2 ) * ( L ( 1 » 2 * 1 i>Ni''l2 ) -5«-U ( 1 » 1 iNM2 ) 

+L t 1 » 2 »2 *imM2 ) I 1 * 2^iMM2 ) +H 1 » 2 » 3i>NM2 ) *U ( 1 » 3it^NiVj2 ) +u I 1 » 2 ♦ 4^Ntvi2 ) 

‘i U ( 1 » 4^i'Ji'^2 ) ) ^ 

up ( 2 * I ^iNii ^2 ) =uP { 2 » I i>NM 2 ) - w*iv, I ( 2 » I « 2'^Nl''i2 ) * < N ( 3 » 2 ♦ 1 ipNKi2 ) *U < 3 ’ 1 it.Nivi 2 ) 

'i +N ^3*2’ 2^i'iM2 ) ^ u I 3 * 2^i'‘i''i2 ) +i'^ ^ 3 * 2 ’ 3^iNi''i2 ) *u C 3 ♦ 3^*'^I''I2 ) +i'J ^3*2* 4^i'^i''i2 ) * 

U ( 3 t 4^Ni^^2 ) ) 

UP ( 2 » I iN,^12 ) = JP ( 2 » I ^NM2 ) - I C 2 » 1 « 3^NM2 ) "^ I L ( 1 » 3 » 1 it»NP ,2 ) *U ( 1 * 1 iNiW2 ) 

^ +L ( 1 ’ 3 » 2^NiVi2 ) "ru ( 1 » 2^i-JPi2 ) +L (1*3’ 3^iMM2 ) *u ( 1 » 3ipN.W2 ) +L (1*3 * 4 ^N.vi 2 ) * 

Lb U ( 1 1 4^imM2 > ) 

UP (2*1 ifclMi'12 ) =UP ( 2 * I -t-N|Vi2 > -W-M-Ivi 1(2*1* 3^1'4M2 ) * ( N ( 3 * 3 * 1 ^^^12 ) "^U ( 3 * 1 i>NM2 ) 

i +N ( 3 * 3 * 2if Nm2 ) *u I 3 * 2i.iNM2 > +N (3*3* 3-J>INi''i2 ) *U I 3 » 3ibiNivi2 ) +.M ( 3» 3 * 4^NiV|2 ) * 

S U ( 3 » 4itNM2 ) ) 

UP (2*1 -liiM'^12 ) = UP ( 2 * I ii^NKi 2 ) 1(2*1* 4^NPi2 ) * ( L ( 1 » 4 * 1 *Ni '^2 ) *U ( 1 * 1 ii;N(''i 2 ) 

i> +L (1*4* 2i-MM2 ) -;(-u t 1 * 2^I'^Pi 2 > +1- (1*4* 3^Ni'^i2 ) *U ( 1 ^ 3ii^NPi2 ) +l_ ( 1 * 4 * 4ii''Mi''i2 ) * 

4 U( 1 . 44NM2 ) ) 

UP (2*1 4NM2 ) =UP ( 2 » I ^Nivi2 ) - I ( 2 * I * 4^iNfvi2 ) * ( iM ( 3 * 4 * 1 ^i’>li''i2 ) ^'U ( 3 * 1 i>Niv'i2 ) 

4 +N ( 3 * 4 ’ 2^I'JM2 ) *u ( 3 * 24Is|vi2 ) +N ( 3 * 4 * 34Ni''l2 ) ( 3 * 3-t'f'''*''i2 ) +i'-i (3*4* 44 Imi''i2 ) * 

i U ( 3 » 44iNli'^i2 ) ) 

UP (2*1 4i'^i'i2 ) =UP ( 2 » I 4i\iiV|2 ) H-iV* ( ivi 1 ( 2 * I * 14NM2 ) *U (2*1 4Nr'i2 ) +l'^i 1(2*1* 2^l'^*''i2 ) 
4 ( 2 * 2 -oPi'"i 2 ) +i-i 1 ( 2*1 * 34N.^12 ) ( 2 * 34Nivi2 ) +i/i 1 ( 2 * 1 * 44Ni'^2 ) *U ( 2 * 4^Nt''t2 ) ) 

UP ( NODE-J » I ) = ( 1 • - vV ) *U ( NODEP * I ) 

uP ( NODEu » 1 )=UP(iV(0OE3» 1 )— w*i^il ( iNODEP * I * 1 ) * ( L ( NM 1 »i»i)*U(Ni-'ll»l) 

if, +L ( NM 1 * 1 « 2 ) U ( 1 * 2 ) +E ( NM 1 ♦ 1 » 3 ) *U ( NiVi 1 » 3 ) +L ( NM 1 » 1 » 4 ) 

4 *U ( NM 1 * 4 ) ) 

UP ( NQDE3 ♦ I ) =UP ( iMODE4 « I ) -iw*M I ( isiUDE3 * I • 2 ) * ( E ( Ni'^i 1 *2*1) *U ( NM 1*1) 

4 +L ( NM 1 * 2 « 2 ) * E ( NM 1*2) +L ( NM 1 *2*3) *U ( Ni^i 1 » 3 ) +E ( NM 1 *2*4) 

4 -;^-u ( NM 1*4)) 

UP ( InODE::^ » I ) = Up ( imUUEE 1 I ) -w*M I ( NUDE3 * i » 3 ) * ( E ( Ni»i 1 * 3* 1 ) *U ( Niv| 1 » 1 ) 

4 +t_ ( Nim 1 * 3 » 2 ) * '“J ( Nim 1*2) +*— ( NMi 1 *3*3)*U( Ni''1 1*3) +E ( NMi 1 *3*4) 

4 *U ( NM 1 * 4 ) ) 

UP ( NOUC-U » I ) = up ( imUUEP * I ) - 1 ( NUUE3 ♦ 1 » 4 ) * ( E ( Ni''l 1 » 4 » 1 ) *U ( NM 1 « 1 ) 

4 +L ( Ni^l 1 ♦ 4 , 2 ) * u ( IMI'^I 1 » 2 ) +E ( Ni'/| 1 * 4 * 3 ) *u ( NM 1*3) +E ( NM 1 » 4 ♦ 4 ) 
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Table D-2 (Concluded) 


* *U ( Niv] 1 » 4 ) ) 

vjP(NODEij» I ) =UP( iMODES* I )+i-v*(|VlI < NODE^ » I * 1 )*D( NODES* 1 )+MI (NODES* I *2 
it ^D ^ NODES* 2 ) +i^ll I i\iOD£S * I * 3) *D ( NODES * 3 > +I''U ( NODEO *1*4) *D( NODES* 4 ) > 

100 continue 

DO 350 1=1*4 

DIF( 1 * ISNODeS) = UP( 1 * ISNODES)-U( 1 * I ANODES ) 

350 continue 
PMS=0 • 

DEL=0 • 

DO 360 1=1*4 

DO 360 J= I * NOoES 

DEL = DEL+D I F ( J > I ) -«-D I F ( J * I ) 

pivis = R|ViS+UP ( J * i ) <^UP( J* I ) 

360 continue 

del=sqpt.( del ) 

RNs=:SOPT (RiViS ) 

TEST=DEL/RNs 
WRI TE ( 6*450 ) UP > J *D IF 
450 format (UE l 0,3 } 

WRITE(6*500) 1 TER*TEST*D£L*RMS 

500 F0RMAT( 15* 3E10 . 3 ) 

DO 400 1=1*4 

U( 1 * I ANODES) =UP ( 1 * ILNODES ) 

4u0 continue 

IF ( TEST.LE.eP5 .OK. I TER. GE. MAX I ) RETURN 

GO TO 10 

END 
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Table D-3 

LISTING OF SUBROUTINE EQSOL5 


subroutine eUSOl!? (U»L *IVUN*D« nodes* w ♦ EPS* MAX I ) 

DIMENSION U(9*5) >L(9*b*5) *M(9«b*5) »MI (9«5»5) * 

N(9’5»5) *DC9*5) »DET(9) » Up ( 9 * 5 ) *DlFt9*5) 
read L*M*MI,|N1 

common/xmi /M l 

M I ( 1 * 1 * 1 a>NODE-^ ) = I M ( 1 * 2 * 2 ^I'^UD£S ) *i''i ( 1 * 3 » 3 SNUDES ) -M t 1*3* E^jIMUDES ) * 

M ( I * 2 « 3-fcNODES ) ) * ( M ( 1 , 4 * 4SNODE^ ) *i''l < 1 « 5 ' 5i*»NODES ) -M (1*5* 44 NUDES ) * 
* M( 1 *4*5SN0dES ) ) 

MI ( 1 * 1 * 1 anodes ) = MI ( 1 * 1 * ISNODES ) - ( M ( 1 * 2 * £ii>NODES ) -M-M ( 1 *4*3SNODES) 

S -M (1*4’ 2 *N 0 D£S ) *i'A ( 1 * 2 * 3 iN 0 DES ) ) * ( M ( 1 * 3 ’ 4 ^N 0 DES ) * 1-1 ( 1 *S * 5 ^faNODES ) 

S — M (1*5* 4SNU0ES ) *M (1*3’ 5 S>NODES ) ) 

Mi ( 1 » 1 * liNOoES) =m 1 (1*1* lSNODLS) + (M( 1 * 2 ’ E^NODLo ) *M ( 1 *5*3SNOUE^) 

S> -M ( 1 * 5 ’ 2S>iNlQL)LS ) *M (1*2’ 3*NODES » ) * ( M (1*3’ 4it»NODES ) *M ( 1 * 4 * 5i»NOUES ) 
S -M ( 1 * 4 * 4iNODES ) *M ( 1 * 3* 5SNUDES ) ) 

M I ( 1 « 1 * 1 -pNud£C 3 ) = 1 - 1 1 ( 1 * 1 * 1 4 jNUDE'^ > + ( I'l ( 1 ’ 3 ’ 2 i^NUDE‘:> ) *M ( 1 * 4 * 3 SNODEO ) 
io “M ( 1 ’ 4 ’ 2^:'N0iJIl.S ) *M (1*3’ 3iNUDE'^ ) ) * ( M ( 1*2’ 4^NODEO ) *M (1*5’ S^NOUES ) 
i -M( 1 *5*4iNODES) *M( 1*2* 5SN0DES) ) 

M I ( 1 * 1 * 1 iNOoEs ) = 1 ^ 1 1 ( 1 * 1 * 1 SNODES ) - ( ivi ( 1 * 3 * 2 ^N 0 DE 0 ) *M (1*5* 3 it>NUDES ) 

^ -M ( 1 ’ 5’ E^NODES ) *M (1*3* 3i&NODE^ ) ) * ( ivi ( 1 * 2 * 4 SNUPES ) *M ( 1 * 4 * SibNODES ) 
it “M( 1 *4*4StNO0ES) *M( 1 *2*53>N0DES) ) 

M I ( 1 * 1 * 1 iNODES ) =i-U ( 1 * 1 * 1 ibNODES ) + ( M ( 1 * 4 * EiJ^NODEi^ ) ■Jti'/i ( 1 * 5 * 3 iti\iUDES ) 
it “M ( 1 ’ 5*2itN0DES ) *M ( 1 * 4 * 35>NODES ) ) * ( M ( 1 * 2 ’ASNODEO ) *M ( 1 * 3 ’ 53.NODES ) 
S “M( 1 *3*4SN0D£S ) *M( 1*2* SitNODES ) ) 

M I ( 1 *2 * 1 itNOpES ) = ( ivl ( 1 * 2 * 1 itNODES ) *M (1*3* 3itNUDES ) -M ( 1*3*1 itNODES ) * 
it> i-it 1 *2*3itiNUDt£C>) ) »((M( 1 *4* 4 i°NODEC :3 ) *l''i (1*5’ Si^NOpto ) ~<*l ( 1 * 5 * 4 5£>N UPE ) * 
it M ( 1 , 4 *5SNOdES ) ) 

M I ( 1 * 2 * 1 itNODE'ii ) I ( 1 * 1 * 2 iNODES ) - ( M ( 1 * 2 ’ 1 itNODES ) *i'/i ( 1 * 4 * 33iN0DE^ ) 

S -M ( 1 » 4 * 1 itNODES ) *M ( 1 * 2 * 3itNOPES ) ) * ( fvi ( 1*3’ 4itNODE:3) *M (1*5* S^NOPES ) 
i -M( 1 ’5’4$N0D£S) *M( 1 *3* SitNODES) ) 

Mi { 1 *2 * 1 SNOpES ) =.^ I ( 1 * 1 * 2itNODES > + ( M ( 1 * 2 ’ 1 i»NODES ) *M ( 1 * 5* 3itN0DES ) 
it -M( 1 ’ 5 ’ litNoDuS ) -i<-M ( 1*2’ 3itN0DES ) ) * (M ( 1 * 3 ’ 4±>NOPEi5 ) *M (1*4’ 5i»NOpES ) 
S “M( 1 *4*4 SNoDES) *M( 1 *3*5SNODES) ) 

MI ( 1*2’ litNOpES ) = Ml ( 1 * 1 * 2itNODES ) + < M ( 1 *3’ l-iNODES ) *M ( 1 * 4 * 3 SN 0 DES ) 

S -M ( 1 ’ 4 ’ 1 itNODES') ^<-M ( 1 * 3 ’ 3itNODES ) ) * ( M ( 1 * 2 ’ AitNOPES ) ( 1 * 5 * SipNODE^ ) 

S ~M( 1 *5*4SN0DES) *M( 1 * 2 * SitNODES) ) 

MI ( 1 *2 ’ 1 SNOpES ) = MI ( 1 * 1 * 2itNODES ) - ( M( 1 * 3 * 1 itNODEi^ ) *M ( 1 * 5 ’3SNODEi:3 ) 
it -’M( 1 ’5’ 1 itNODES) *M( 1 *3’3SN0PES) )*(M( 1 * 2 ’ 45.NOPES ) *M ( 1 * 4 ’SitNODES) 
it -M ( 1 * 4 *49 MODES ) *M ( 1 * 2* SitNODES ) ) 

M I ( 1 * 2 ’ 1 9 >NOpES ) = M I ( 1 ’ 1 ’ 2 itNODEii> ) + ( M ( 1 * 4 ’ 1 ibNODEis ) *i''i ( 1 ’ 5 ’ 39 >NOpEi 3 ) 
it -M ( 1*5*1 ioNODEO ) (1*4* SitNODES ) ) * ( M ( 1 * 2 • 4ibNOPES ) *M (1*3* SitNUpEi^ ) 

S. -M( 1 *3’4^N0DES) *M( 1 * 2 * SitNODES) ) 

M I ( 1 * 3 ’ 1 itNOpES ) = ( M ( 1 * 2’ 1 itNODES ) *M (1*3’ 2itNupES ) -M ( 1 * 3*1 itNODES ) * 

(Continued) 
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Table D-3 (Continued) 


it M ( 1 ♦ 2 » 2itNODE,S ) ) * ( M ( 1 * 4 » 4itN0DLS ) *|vi (1*5* SitNODtii:^ ) ~M (1*5* 4itNUDEi3 ) * 

^ M( 1 * 4*5SNOdES ) ) 

( 1 * 3 * 1 i*>NODE^ ) =i'11 ( 1 * 1 * ) - ( M ( 1 * 2 » 1 itoNODh.i:> ) (1*4 * 2iNuDE3 ) 

it -M( 1 *4» li&i'iODtb) *i'^i( 1 *2»2itN(jDES) ) * ( M ( i * 3 » ^iNODES ) ( i *5*5itNQDES) 

S -M ( 1 * 5» 4itl\l0DES ) (1*3* SiNODEE > ) 

M I { 1*3*1 SNGdE-S ) = 1 '^ I ( 1 * 1 * 3it>N0DES ) + (M( 1*2*1 SNODEb ) *(V) ( 1 * 5 * 2$NODES ) 
it -jvi ( 1 , 5 , ISNQDES ) ( 1 * 2 * 2 itNOD£S ) ) * ( M ( 1 * 3 * 4±N0UES ) *M ( 1 * 4 * SiNODES ) 

it -M (1*4 *4itNQDEb ) (1*3* 5S>N0DEb ) ) 

I ( 1*3*1 it»NOoE-ii> ) =i-'l 1 ( 1 * 1 * 3 it>NuuE:D ) + ( I'^i ( 1 * 3 * 1 :ti\ODE-> ) *i'/l ( 1 * 4 * 2 itiNuD£b ) 

it — |V| ( 1*4*1 *NODt_C:3 ) ( 1 * 3 * 2^NODEb ) ) * ( iv'i ( 1 * 2 * 4it>N0DE3 ) *M (1*5* SitNijDEb ) 

it -,'Vl ( 1 * 5* 4iN0DcS ) ( 1 * 2* SibNODES ) ) 

M I ( 1*3*1 itNOoEb ) =i'^l I ( 1 * 1 * 3it>NODEb ) — ( ivi ( 1 * 3 * 1 ibNUDE-^s ) *l '>1 ( 1 * 5 * 23 >NODEi> ) 
it -|V)( 1 * 5 * lii3iN0Diio) *J''i( 1 * 3* 2it'N'JDEE ) )*(i'vl( 1 *2 »4itN0DEb) ( 1 *4*5iNUDEb) 

i -M( 1 *4 *4ii>N0DES J 1 * 2 * 5S:>NODEb > ) 

M I ( 1 * 3 * 1 itNOD ) = i''tl ( 1 * 1 * 3 itNODEG ) -♦- ( M ( 1 * 4 * 1 iNODE^ ) ( 1 * 5 * 2itN0DEii^ ) 

it ~M ( 1*5*1 itNODlib ) ( 1 * 4 * 2 itNOD 5 b ) ) * ( |V) ( 1 * 2 * 4^NUDEb > *M (1*3* S^tNODES ) 

i -M( 1 *3*4itN0DES ) *^( 1 *2* SitNODES) ) 

M 1 ( 1 * 4 * 1 ibNUDEii> ) = ( M ( 1 * 2 * 1 itNODEb) ) (1*3* 2itNOUEb ) -M ( 1*3*1 itNUDEb ) * 

» .vl ( 1 * 2 * 2ibNUDE^ ) ) ^ ( I''! (1*4* 3iwiNGUEb ) *|V| (1*5* 5ibNUDE^ ) ( 1 * 5* 3itNODE^ ) * 


a. i'i( 1 *4*5iiNODEb ) ) 

M I ( 1 * 4 * 1 itNUDEiiJ ) =M I ( 1 * 
a. -|Vi( 1 * 4 * litNODEb ) ( 1 * 

S -M ( 1 * 5* 32>N0DES ) ( 1 * 

IV, 1 ( 1 , 4 , 1 4>nud£^ ) =M I ( 1 * 
it -R ( 1*5*1 itiNQDES ) ( 1 * 

ib -R ( 1 * 4 * 3itNODES ) *R ( 1 * 

R I ( 1 1 4 * 1 3d E M I ( 1 * 

i — R ( 1*4*1 it.NiODEb ) ( 1 * 

it -M ( 1 * 5*3itN0DES ) ( 1 * 

R I ( 1 * 4 * 1 ipNODE-^ ) =i"'i 1(1* 
it -R ( 1 * 5* 1 itNQOEb 1 * 1^1 ( 1 * 
i -R ( 1 *4*3a>N0Dtb ) -i«-M( i * 
Ri ( 1*4* l^NODEb) =MI ( 1 * 
it ~R( I *5* I it NOD Eb ) *R( 1 * 
i -M ( 1 * 3 ♦ 3itNODES ) *R ( 1 * 
R I ( 1 * 5 * 1 itNODEb ) = ( M ( 1 * 

St R ( 1 * 2 * 2itNCJD£o ) ) * ( M ( 1 
2> M( 1 * 4*4£NOdES ) ) 


1 * 4itNODEb ) - ( R ( 1 * 2 * 1 itNODEb ) *R (1*4* ZiNUDEb ) 
2* 2-'tN00Eb) )* (R( 1 » 3 * 3 itNODEb ) *M ( 1 * 5 * 5 *NODEc 3 ) 
3* SitNODES ) ) 

1 * 4ibNODEb ) + ( R ( 1 * 2 * 1 itNUDEiib ) *R ( 1 * 5 * 2 itNQDEb ) 
2* 2itNUDEb ))■»■( R ( 1 * 3 * 3^NUDEb ) *R (1*4* 5*N0DEb ) 
3* SitNODEb ) ) 

1 * 4a>N0D£b ) + ( R ( 1 *3 « 1 SNODEb ) *M ( 1 * 4 * 2 itNODEb ) 
3 * 2ibNUDEi::5 ) ) * ( R (1*2* 3 itNODEb) *M (1*5* SibNEDEb ) 
2 * SSNODEb ) ) 

1 * 4i»N0DEG ) - ( R ( 1 * 3 * 1 ibNODi^'b ) *R ( 1 * 5 * 2itNUDEi=> ) 

3 * 2 itN 0 DEb ) ) * ( R ( 1 * 2 * bibNUUEb ) *R ( 1 * 4 * SiNODES ) 

2 * SibNODEb i ) 

1 * 4ibNUDEb ) + ( R ( 1 * 4 * 1 itNUDEb ) *M ( 1 * 5 * E^NODEb ) 

4 * 2 itNODES ) ) * ( M ( 1 * 2 * 3 itN 0 DEb ) *M ( 1 * 3 * SitNUDEb ) 

2 * SitNODES) ) 

2 * 1 itNODEb ) *M (1*3* 2itN0DEb ) -R ( 1*3*1 itNODES ) * 

* 4 * 3SbNODE^ ) *R (1*5* 4 itNUDEb ) -R (1*5* 3 itNOUEb ) * 


R I ( 1 ’ 5 * 1 itNODES ) = 1 ^ 1 1 ( 1 * 1 * 5 i»NODEb ) - ( R ( 1 * 2 * 1 itNODEb ) *M ( 1 * 4 * 22 .N 0 DES ) 
it -R ( 1*4*1 itNODEb ) ■’'^R (1*2* 2 ibN(JDEb) ) * ( R ( 1 * 3 ' 3 itNODEb > *M( 1 *5 » 4 itNODEb ) 
it -*M( 1 *5*3iN0DES } ^M( 1 *3* 4itN00ES) ) 


R I ( 1*5*1 itNODEb ) =i'l 1(1*1* SitNODEb ) + ( R ( 1 * 2 * 1 i»NUDE^ ) *R (1*5 * 2itNOOEib ) 
ib “R (1*5*1 ibRQDEs ) ( 1 * 2 * 2^>'^0DES ) )*(r( 1*3* 3itROOEi:> ) *i''i (1*4* 4ibi'40DEb ) 

it “R ( 1 * 4 * 3it NODES ) *R ( 1 * 3* 4itNODEb ) ) 

R I ( 1 * 5 * 1 ibNODEs ) =t"U ( 1 * 1 * SibNUDEb ) + ( R ( 1 * 3 * 1 ibNuDES ) *R (1*4 * 2 itNUDEib ) 
ib ~R( 1 * 4 * libNODi^S) *R( 1 * 3 * 2 itNODEb) ) -)t ( R ( 1 » 2 * 3itN0DEE ) -«rR ( 1 *S*4ibNODEb) 
i ~M( 1 »5*3 SNoDES) *R( 1 *2*4SN0DEb) ) 

(Continued) 
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Table D-3 (Continued) 


i''il ( 1 «5* li*iNUDES) = ,vll { 1 * 1 * 5 *N 00 ES)“ ( M( 1 1 ibNODEiJ ) ( 1 » 5 * 2 itNODES) 

3t> t 1 ’ 5* liNODEia ) <1*3* 2 ^NOQEb ) ) * ( iVi ( i , 2 *3^NODtij ) *M ( 1 « 4 « 43>NODE3 ) 

S -M< 1 »4*3^N0DES ) *M( 1 *2*4aNODES) ) 

MI < I *5* ISNOdES J =MI ( 1 • 1 • 5^£>N0DES) + (M( 1 *4» 1 SNODEb ) *|V| ( i » 5 * 2 ^NODES) 

S -Mt 1 *5’ l^tNODbS) *M( 1 *4*23>NODEi>) ) * ( M ( 1 » 2 ’ 3*N0DE3 > *M I 1 <3»4*N0DES) 

S “ivi < I ♦ 3 • 3SNODES ) *M ( 1 « 2 » 42>N0DES ) ) 

M I < 1 * 1 » 2^i»N0dE^ ) = ( M ( 1 , 1 , 23tNODE^ ) *M ( 1*3* 3^tN0DE^ ) ~M ( 1*3* 24>NuDE3 ) * 

S> M( 1 * 1 *3^NOdE 3 ) ) * ( M ( 1 *4 * 4ibNODES ) *M( 1 * 5 * S^NUDES ) -M ( i * 5, 4i&NQDES ) * 
M( 1 *4*5*NOdE3 ) ) 

MI ( 1 * 1 *2$N0dE^)=MI (1*2* ia>NODES)- { M( 1 * 1 * 2SNODLS ) *|Vl ( 1 *4*3SN0DE3) 

S. <1*4* 2 *N 0 DES ) -«-M ( 1 * 1 * 3 i»»NODES ) ) * ( M ( 1 * 3 * 4 ^N 0 DE 3 ) -Jt-M ( 1 * 5 * SibNODES ) 

S -M( 1 *5*4SN0DES) *M( 1 *3*5itN0D£3) ) 

MI ( 1 , 1 *2itNODE3 ) = M I ( 1 * 2* l^NUOE^ ) + C M ( 1 * 1 * 2iNODE3 ) *|>/1 ( 1 * 5 * 3 SNODES ) 

S “M ( 1 * 5* E'i’NODE^ ) *M ( 1 * 1 * 3^NodE3 ) ) * ( M ( 1 * 3 * 4^InUDE3 ) *M( 1 * 4 * 5s»NudE3 ) 

^ -M( 1 *4*43.N0DES ) *M( 1 *3*53>N0DE3) ) 

M I ( 1 * 1 * 2itoNODE3 ) = M I ( 1*2* 1 *N0DE3 ) + ( M ( 1 * 3 * 2 ^±'MODE 3 ) *)Vi ( 1 * 4 * 3 SNQDE 3 ) 
it. -iV] C 1 * 4 * 2 ^NODii 3 ) *M ( 1 * 3* 3^N(JDE^ ) ) * ( M ( 1 * 1 * 4 *N 0 DES) *M ( 1 * 5 * 5 ^*»N 0 DES ) 
it -M ( 1 * 5*4iBN0DES ) *M ( 1 * 1 * 53>N0DES ) ) 

MI { 1*1 ♦2it^ JdE^ ) = .^l I ( 1 * 2* litNODE^ ) - ( M ( 1 * 3 » 2it>NUDEi> ) *M ( 1 ♦ 5* 3*NODE3 ) 

S “M < 1 *5*2iiM0D!::3 ) *M ( 1 * 3* 3iNODE3 ) ) * ( M ( 1 * 1 *4itNOUE^ ) *M( 1*4* 5SbN0DE3 ) 
it -|Vi( 1 *4 *4itNODES ) *M( 1 * 1 fSitNODES) ) 

M I ( 1 * 1 * 2itNODE3 ) =M I ( 1*2* li&NUDEb ) + ( M ( 1 * 4 * 2itNODE3 ) *M ( i * 5 * 3^NODES ) 
it ~M ( 1 * 5 *2i»MODE3 ) *M I 1 * 4 * 3i»NODE^ ) ) * ( M ( 1 * 1 » 4it'tNKJDE3 ) *M( 1 * 3 * S^NUDEb ) 

S -M( 1 *3*43.N0DES) *M( 1 , 1 ,55.NUDEb) ) 

MI ( 1 *2*2itNUDES) = (M( 1 * 1 * 1 ^NUDES)*m( 1 * 3 * 3itNODEb J -M ( 1 *3* 1 iNODES ) * 
it M{ 1 * 1 *3^N0dES ) ) *( M( 1 , 4 * 4itNUDES) *M ( 1 *5*5i»NODEb ) -M ( 1 * 5 » 43 >NUDEb ) * 

S M{ 1 *4*5itN0DES ) ) 

M I ( I * 2 *2itN0DE3 } = M I ( 1 • 2* 2itNODES ) - C M ( 1*1*1 SNODES ) *M ( 1 *4 * 3 ^N 0 DES ) 
it -M( 1 , 4 * 1SN03E3)*M( 1 , 1 *3itN0DES) )*(M( 1 *3*4iNUDES) *M ( 1 *5*5itN0DES) 
it ~M( 1 * 5»4itN0DES ) *M( 1 *3* SitNODEb) ) 

M 1 ( 1 * 2 * 2itNooEb ) 1 ( 1 * 2 * 2^MUDEb ) 4- ( M ( 1 * 1 « 1 itNODEb ) *M (1*5* 3itNUDEb ) 

it — m( 1*5*1^ NODLb) 1 1 , 3itoNUDEb ) ) * ( ( 1 * 3 » 4i>NODE':s ) *M ( 1 *4*5iaNUDEb) 

i -M( 1 *4»43 sNODES ) +^M( 1 * 3* SitNODES) ) 

MI ( 1 *2 » 2 itNODEb ) = M I ( 1 * 2* 2 itNODEb ) + ( M( 1 *3* 1 itNUDEb ) *M ( 1 * 4 * 3 ^NODEb ) 
it ~M ( 1 *4 ♦ 1 itNOOtb ) (1*3* 3iNOQEb ) ) * ( M ( 1 * 1 * 4ifaiMUDEb ) *M ( 1*5* SitMOOEb ) 

S> -M ( 1 *5»4itNODES ) *M (1*1* 55.NODES ) ) 

M I ( 1 * 2 » 2SNOdE^ ) =M I ( 1*2* 25>NODEb ) - ( M ( 1 * 3 * 1 itNODEb ) *M ( 1 * 5 * 3 itNODEb ) 
it *-M ( 1*5*1 SNODE3 ) -Jt-M ( 1 * 3 * 3ibN0DE3 ) ) * ( M ( 1 * 1 * 4S>N0DEb) *M ( 1 * 4 * SitNODEb ) 

S -M( 1 *4*4SN0DES ) -«-M( 1 * 1 * SitNODEb) ) 

MI ( 1 *2*23>NQdES) =MI ( 1 *2* 2 itNODEb ) + ( M ( 1*4* 1 SNODEo ) *M ( 1 *5*34>N0DEb) 
it — M ( 1*5*1 itNODEb ) *iM ( 1 * 4 * 3 ifNODEb ) ) * ( M ( 1 * 1 * 4^NODEb ) *M (1*3* 5 i*»NUDEb ) 
i ~M( 1 *3*4SN0DES ) *M( 1 * 1 iSSNODES) ) 

MI ( 1 *3*2SN0dE^ ) = (M( 1 ♦ 1 * 1 itNODEb ) *M( 1 * 3 * 2i^NODEb ) -M ( 1*3* liNODEb )* 
it M (1*1* 2it>NODE'^ ) ) * ( M ( 1 * 4 * AitNODEb > *)vi <1*5* SibNUDEb ) -M (1*5* 4 *NODEb ) * 
it M ( 1 « 4 * SitNODEb } ) 

MI ( 1 *3*2SNOdES) =MI ( 1 *2*3itN0DEb)-(M( 1*1* 1 SNODEb ) *M ( 1 *4*2SN0DEb) 
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Table D-3 (Continued) 


5 b -iVi C 1 *4 » 1 anodes ) ( 1 * 1 » SiNOOEia ) ) *' ( M ( 1 *3 *4*NODES ) ( 1*5* SiNODES ) 

ib “M( 1 * 5 » 43 >N 0 DES) 1 tafS^bNODES) ) 

MI ( 1 *3’2ibNQDEo ) =MI ( 1 1 2 < 33>NODEb ) + ( M t 1^1* 1 5bNUDE^ ) ( 1 « 5 * aa^NODEia > 

S -Mi 1 »5* 1^1 'JODE 6 ) *M ( 1*1 »2a>NODES) )*(M( 1 * 3 * 45 bNODEi ) ( 1 *4*5iNUDEii> 

S -M( 1 »4 *4a:;MODtiS ) ( 1 *3*5^N0DEE) ) 

MI ( 1 *3«2iN0DE^) =MI I 1 « 2 * 3 S>NODES) + (iVU 1*3* 1 SNODES ) -k-M ( 1 »4»2iN0DES) 

S -M ( 1 *4 * ISNoDtS ) *M ( 1*3* 2 SNODES) ) * ( M ( 1 * 1 *4SN0DES) *M ( 1 *5* 53>N0DES ) 

-M ( 1 *5*4SN0DES) -^M( 1 * 1 *5a>NODES) ) 

M I ( 1 * 3 *2SNOdE*:> ) =M I i 1 * 2 * 3 ibN 0 DE^ ) - ( M i 1 *3 * 1 ibNODE^ ) * 1^1 i 1 * 5 * 2 SMODEi:> ) 

4 -Mtl*5*ia*N0Dt3) *M i 1 * 3 * 2a>NODE3 ) ) -K- ( M ( 1 * 1 * 4a5NCDE3 ) -X-M ( 1 * 4 * 5i>NOD£3 ) 

5 -M( 1 *4*4^NODhS) *M( 1 * 1 *53^N00ES) > 

MI ( 1 *3*24iMOd£3) = MI ( 1 * 2* 3^NUDEE ) + ( M ( 1*4* 1 4NODE3 ) ( i *5*25>N0DES) 

i -M ( 1 * 5 * 1 SNODES ) *M (1*4* 23=N0DES ) ) * ( M ( 1 * 1 * 4SNODE3 ) *M ( 1 * 3 * 5*N0DES ) 

S -M< 1 *3*4SN0DES ) *M( 1 * 1 *5SNOOE5) ) 

M I ( 1 * 4*2SN0dEE ) = ( M { 1 * 1 * 1 SNODES ) *M ( 1 * 3 * 2a>N0DES ) -M ( 1 * 3 * 1 SNODES ) * 

S M ( 1 * 1 * 2SNOdE^ ) ) ^ ( M ( 1 *4 * 3 a>NODES ) -M-M (1*5* 5SNODES ) -|v| (1*5 * 34NODES ) * 
S M( 1 *4*5SNOdES ) ) 

M I ( 1 *4 *2SfV0DES ) =M I ( 1 * 2 * 4SNODES ) - ( M ( 1 * 1 * 1 SNODES ) *M (1*4* 2^bN0DES ) 

S -M( 1 *4* ISNODES) *M ( 1 * 1 *2^N0DES) )*(M( 1 * 3 * 3SNODES ) *M ( 1 *5*5SNODES) 
S -M( 1 *5*3SNoDES) *M( 1 *3*5SN0DES) ) 

MI ( 1 * 4*2S,N0dES ) =MI ( 1 * 2 * 4SN0DES ) + ( ( 1 * 1 * 1 SNODES ) *M ( 1 * 5 * Z^NODES ) 

S -M (1*5* ISNODcS ) *M ( 1 * 1 * 2 SNODEi 2 ) ) * ( IM ( l ♦ 3 * 3 SNOCES ) *M ( 1 * 4 * S^^NODES ) 

S -M( 1 *4*3SN0DES ) *M( 1 *3*5SN0DES) ) 

M I ( 1 * 4 *24N0dES ) =i-4 I ( 1 * 2 * 4 SNODES ) + ( M ( 1 * 3 * 1 SNODES ) *M (1*4* 24N0DES ) 

4 -M ( 1 * 4 * 1 SNCDES ) *M ( 1 * 3 * 2^N0DES ) ) * ( N ( 1 * 1 * 3 SNODES ) *M ( 1 * b * 54NODEO ) 

5 “*M( 1 *5*3SNoDES ) *M( 1 * 1 *5SN0DES) ) 

MI ( 1 *4*24N0dES) =MI ( 1 *2*44NoDES)~(M( 1*3* 1 SNODES ) *M ( 1 *5*23>N0DES) 

S -M ( 1 * 5* 1 SNODES ) -i<-M ( 1 * 3* 2 SN 0 DES ) ) * ( M ( 1 * 1 * 3 SNODES ) *M ( 1 * 4 * SSNODES ) 
S -M( 1 *4*3SN0DES) -«-M( 1 * 1 ,5SN0DES) ) 

MI ( 1 »4*2SN0dES) =MI ( 1 » 2 « 4 SNODES) + (|Vi( 1 * 4 * 1 SNODEi^ ) *M ( 1 *5*2^bN0DES) 

S -M ( 1*5*1 SNODEE ) -«-M ( 1 * 4 * 2 SNODES ) ) * ( M ( 1 * 1 * 3 SNODES) *M ( 1 * 3 * 5SN0DES ) 
S -M ( 1 *3*3SN0DtS ) *M ( 1 * 1 * 5SNODES) ) 

M I ( 1 * 5* 2*N0dEo ) = ( M ( 1 * 1 * 1 5bNODES ) *M (1*3* 2^NODE-:> ) ~M ( 1*3*1 SNoDES ) * 

S M ( 1 * 1 * 2SN0dES ) ) ^ ( M ( 1 * 4 * 3 SNODES ) (1*5* 4 SNODES ) -M ( 1 * 5 * 3 SNODES ) * 

4 M ( 1 * 4 * 44NODE5 ) ) 

MI ( 1 »5*24NOdEo) =1^11 ( 1 * 2 * 5SN0DES ) - ( M ( 1*1* 1 SnODE:? ) *1^1 ( 1 *4*2SN0DES) 

5 -M( 1 *4* 1 SNODES ) *Mi( 1 * 1 * 2 SNODES) >*(M( 1 * 3 * 3SNQDES ) *M ( 1 *5*4SN0DES) 

S -M ( 1 * 5* 3 SNODES ) (1*3* 4SN0DES ) ) 

MI ( 1 *5*24NOd£S) = MI ( 1 * 2 * 5 SNODES) + (M( 1 * 1 * 1 4NODE3 ) "'(‘M ( 1 *5*2SN0DES) 

S -M ( 1 * 5* 1 SNODES ) ( 1 * 1 « 2 SNCDES ) ) * ( M ( 1 * 3 * 3SNODES ) *M ( 2*4 t 4 SNODES ) 

4 -M( 1 *4*3SN0DES ) -«-M( 1 *3*4SN0DES> ) 

I ( 1 * 5 *2SN0dES ) =i 'l 1(1*2* S^NODES ) + ( jvi ( l * 3 * l 4NODE3 ) *M (1*4 * 2 SNUDES ) 

4 — M ( 1*4*1 anodes ) (1*3* 2^M0DE*^ ) > * ( i''i ( 1 ♦ 1 * 34(NiUDES ) (1*5* 4-<^NuDES ) 

4 -M ( 1 * 5 * 3SNODES ) *M (1*1* 4 SNODES ) ) 

MI ( 1 *5*2^i'J0DE'^) ( 1 *2* SSinODE^ )~(M( 1 *3* IsNODE^) (1*5* 2 SNODE^ ) 
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* 1 «5* lS.NOUtS) *M( 1 <3»2S>NQDES) 1 » 1 O^NODES)*MC 1 »4i4^NC0£S) 

S -M ( I * 4 ♦ 3*NoDES ) *M ( 1 * 1 » 4a>N0DES ) ) 

Ml ( I 2i.NOD£i> ) = MI I 1 * 2* S^i'NOOES ) + ( M( 1 -i 4 » 1 ii.NODEb ) -5t-M ( I »5*2^fcNODES) 

^ -M( 1 «5» I^^NODlS ) *M( 1 * 4 * 2^N0DEi=> ) )*(M< 1 » 1 «3^N0UES)*M( 1 ♦ 3 * 4i»>N0DES ) 
^ “M( 1 ♦3^3iNoD£S 1 *M( 1 , 1 <4^fcNOOES) ) 

MI ( 1 » 1 « 32-NOdES ) = ( M ( 1 » 1 « 2^NQDES ) *l^i (1*2* 3 ±NOdE£> ) ”M( 1 * 2 * 23>N0DES ) * 

S M( 1 * 1 *3ii>NODES) )*(M( 1 * 4 * 4ifcNODES) *M( 1 . 5 * S^tNUDES ) ~M ( 1 » 5 * 4 SNODES ) * 
S M( 1 »4*5£N0dES ) ) 

MI ( 1^1 *3SN0dEC)) =MI (1*3* 1SN0DES)~(M( 1 * 1 » 2 *NODEc 5 ) *M ( 1 * 4 * 35 .N 0 DES ) 
ifc -|V 1 ( 1 ♦ 4 * 2 ^N 0 DES ) ( 1 , 1 , 34 >NQL)Ei 3 ) > * ( M ( 1 * 2 * 4 ^MUDES ) ( 1 * b * 5 S>N 0 DES ) 

S -M( 1 *5*4ai\oDE5 ) *M ( 1 ^ 2 * S^NODES) ) 

M I ( 1 * 1 * 3itNQDE-i ) =MI ( 1*3* lipNODEb ) + ( M ( 1 * 1 * 2ipNUDE::^ ) *M ( 1 * 5 * 3 *NUDEb ) 

a -M( 1 *5*2 ^NoD£3 ) ^«-M( 1 t 1 » 3 ^NUDEb) )*(M( 1 * 2 * 4iNQDEb ) ^M ( 1 »4*5^MUDES) 

S “M( 1 *4*4it>N0DES ) *M( 1 *2* 5itNOOES) ) 

MI ( I * I * 3ifcN(JDE3 ) = MI ( 1 * 3* l^iMOUE:^ ) i- ( M ( 1 • 2 * 2 ^MUOEb ) *M ( 1 * 4 • 3ifcNuDE^ ) 
i -M ( 1 , 4 *2^N0DtS ) *M (1*2* 3^N0DEb ) ) * ( M ( 1 * 1 * 4^ttN<ODES ) ( 1 *5 * S^NODEi) ) 

£ “M( 1 *5*4^iNODtS ) -»(-M ( 1 , 1 ,5iNODES) ) 

M I ( 1 f 1 * 3itN0DEb ) =1^11 ( 1*3* l^fcNODES ) -- ( J^l ( 1 * 2 * 2 itNuDEb ) *IVt ( 1 * 5 * 3 ii»NODE o ) 
i “M (1*5* 2itNODE3 ) *M (1*2* 3^MODE^ ) ) * ( M ( 1 * 1 * 4*N0DES ) *M (1*4* S^MODE:^ ) 
S -M( 1 *4*4^N0DES ) *M( 1 * 1 *53’N0DES) ) 

MI ( 1*1 * 3 ±N 0 dES ) = M I ( 1*3* ia.NODES ) + ( M ( 1 * 4 * 2 i£»N 0 DES > *M ( 1 * 5 * 3 itNODEb ) 

^ -*M ( I * 5 ’ 2 ^M 0 DEb ) *1''! ( 1 * 4 * 3 bNUDEb } ) * ( M ( 1 * 1 * 4^NuD£b ) *ivi (1*2* S^NUDE::^ ) 

* -M( I *2*4iN0DE3 ) *M( 1 * 1 *5^NUDES) ) 

M I ( 1 * 2 * bipNOoE^ ) = ( I'i ( 1 * 1 * 1 ipNUDEb ) (1*2* 3ifaNOu)Eb ) -M ( 1*2*1 iNUDEb ) 

a> M ( 1 * 1 * 3SN0DE>i ) ) * ( M ( 1 *4 * 4^NODEb ) *M ( 1 * 3 * SSNUDEb ) -i^i ( 1 * 5 * 4 ^NODEb } * 
S> M( 1 * 4 * 5ii»N0D Eb ) ) 

M 1 ( 1 * 2 * 3*NQDEb ) = M I ( 1 * 3* 2 =i’NUDEb ) - ( M ( 1*1*1 ipNOD£i> ) *N ( 1 * 4 * 3 SNODE 0 ) 
a, -M ( 1 ♦ 4 * libNODt.^ ) *M ( 1 * 1 * 3 ^NUL;)Eb ) > * ( M ( 1 * 2 * 4 ^NQDEb ) *M( 1*5* S^NODE^ ) 
S -M ( I * 5*43 nIODES ) *M (1*2* S^NODEb ) ) 

Mill * 2* 3iN0DEb ) =M I ( 1 * 3* 23=’NUDEb ) -f* ( M ( 1 * 1 * 1 itNODE^ ) *M (1*5* 3iN0DEb ) 

S -M ( 1*5*1 iNODES ) *M ( I * 1 * 3 ^NODEb ))-«■( M ( 1 * 2 * 4*NQDEb ) *M (1*4* S^MODEij ) 
$ -M ( 1 *4 *4SN0DES ) *M ( 1 * 2 * S^NODES ) ) 

MI ( 1 *2*3SNUDEb ) = MI ( I * 3* 23’NODEb ) + ( M( 1*2* 1 JcNOOEb ) *M ( 1 *4*3a>N0DES ) 
a> -M ( 1*4*1 a»i'JODEb ) ^<-M (1*2* 3i°N(JDEb ) ) * ( M ( 1 * 1 * 4i*>iMOUEb ) (1*5* S^NOUEc^ ) 

i -M( 1 *5*4 SuNOD£ 5 ) *M (1*1 *5^tN0DEb) ) 

M I ( 1 * 2 * 33>NUd£o ) =M I ( 1 * 3* 23»MODEb J - ( M ( 1 * 2 * 1 a^NODES ) -«-M ( 1 * 5 * 3 a>NUDEb ) 

^ -M ( 1 * 5* l^NODEb ) *i’^i < 1 * 2 * 3i>NODEb ) ) * ( M ( 1 * 1 * 4^NUDEb ) -i<-M ( 1 *4 * 53>N0DE:3 ) 

* ~M( 1 *4*4*N0DES) *M( 1 * 1 ,53»NODEb) ) 

Mill ♦ 2*3 SNOdE3 ) =MI ( I ♦ 3 * 23>NODEb ) + ( M( 1 * 4 * 1 S-NUDES ) -J(-M ( 1 *5 *3SN0D£b ) 
i -M ( 1 * 5* 1 anodes ) *M ( 1 * 4 * 33 iNODEb ) ) * ( M ( 1 * 1 * 4 iNODEb ) *|V| ( 1*2* 5 *NODEb ) 
S ~M( 1 *2*43iN0DES ) -«-M( 1 * 1 *5*N0DES) > 

Mil 1 * 3*3a.NOoEb ) := (M I 1 * 1 * ISNODEb) *M( 1 * 2 * 2 *NOOEb ) -M ( 1*2*1 SNODES ) * 

2. Ml 1 * 1 * a^iNOoEb ) ) * I M ( 1 *4 * 43 »NODEb ) -K-M ( 1 * 5 * 53 jNODEb ) -M I 1 * 5 * 4 fbNODEb ) -Jf 

S M I 1 * 4 * 52>MOdES ) ) 

Mill * 3* 3SNODE3 ) =M I I 1 * 3* 33>NODEb ) - I M ( 1 * 1 * 1 SNODEb ) *M ( 1 * 4 * 2 SNODES ) 

S -M( 1 *4* ia.NOD£S) *Ml 1 * 1 * 23»NODEb ) )*(M( 1 * 2 * 4^N0DEb ) *M I 1 *5*53>NODEb) 

* -M ( 1 * 5*4iNOD£S ) *Ml 1 * 2* SSNODES ) ) 
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Table D-3 (Continued) 

MI ( 1 »3»3i>N0DEc>) =MI ( 1 » 3» 3it.N0DE3) + ( M( 1 » 1 < 1 -ANODES ) ( 1 t5*2*N0DES) 

i -jVi( 1 litNOD£3 ) 1 « 1 »2^N0DE3) ) * < |vi< i » 2 * 4^iNUDES ) *|V| ( i «4«5^N0DES ) 

S -M( 1 *4*4*IM0D£S ) 1 f2*5^NODES) ) 

Mill « 3* 3ibN0D£3 ) =rl 1 I 1 < 3* 3 SpNOUE 3 ) + ( M ( 1 ♦ 2 * 1 ifcN0DL3 ) *M ( 1 f 4 « 2^M0DE3 ) 

S -MC 1 »4* IS-NODES) 1 *2»2^N0DE3) )*(M( 1 ♦ l * 4H.N0DES ) -X-|vi ( i »5*5±N0DE3) 

S -Ml 1 *b»4^N0D£S) *M( 1 f 1 tSitNODES) ) 

M I ( 1 1 3« 3iNODE3 ) =M I ( 1 » 3i 3^N0DES ) - ( M ( 1 » 2» 1 ibNODE^ ) -i«-M (1*5* 2 S>NuDES ) 

S “M < 1 • 5 * 1 *NoD£3 ) -?{-M ( 1 ♦ 2 * 2^N0dE3 ) ) * ( M ( 1 « i » 4ifNODEb ) ( 1 < 4 * 5i>NUDE3 ) 

S -M ( 1 *4 *4^i^ 0DES ) *tvi ( 1 , 1 , SibNODES ) ) 

Mill » 3» 3iNUDE^ ) = M I ( 1 » 3' 3ibNUDE3 ) + I M ( l »4 » 1 ^INOQE^ ) *M ( 1 » 5» 2^tN0DEiJ ) 

S -M( 1 »5< lipiMOOti-CD ) *M( 1 »4»2^N0DE3J )*(M( 1 » 1 » 4^bNUUEiJ ) ( 1 t 2 « S^^NUiJEh ) 

3 > -M t 1 * 2 ' 4*N0DES ) *M ( 1 » 1 » 5SN0DE3 ) ) 

MI ( 1 » 4 *3^£'^ 3dE^ ) = ( M ( 1 * 1 * 1 ANODES ) (1*2* 23>NUDES ) ~M ( 1*2*1 ANODES ) * 

± M (1*1 *2ii>N0DEc> ) ) * ( M I 1 *4 * 3SNODE3 ) *M ( 1 * 5 * SiijMODE::^ ) -M ( 1 * 5 * aSNuDE::^ ) 

S M( 1 *4*5iNODES) ) 

MI { 1 * 4 * 3i£>NODEo ) =i''|I ( 1 * 3 * 4ibNODES ) - ( M ( 1 * 1 * 1 iNUDE^ ) *M ( 1 *4*2$NODEi>) 

i -M ( 1 * 4 » liNODEC^ ) >M ( 1 * 1 * 2it>NODE3 ) ) * ( M ( 1 ♦ 2 * 3^t>N0DEb ) ( 1 * 5 * S^NODEb ) 

S -M( 1 * 5* 3ii>NoDES ) *M( 1 *2*5iN0DEb) ) 

M 1 ( 1*4* 3SNODES ) =M I ( 1 ♦ 3* 4a>NODEb ) + ( M ( 1 * 1 * 1 it-NODEb ) *M (1*5 * 2a>NODES ) 
i -M( 1 *5* ISNoDES) ^M( 1 * 1 *2^bNODEb) )■;<■( M ( 1 ♦ 2 * 3^i>NaDES ) ( 1 *4*5±NODES) 

S -M ( 1 * 4 * 3SN0DES ) *M ( 1 * 2* 5SNODES ) ) 

M I ( 1*4* 3i»NUD£o ) =r'1 1 ( 1 * 3 * 4ibN0DEb ) + ( M ( 1 * 2 * 1 iNODEb ) ( 1*4 * 2 *NUDEb ) 

it> ~M ( 1*4*1 i>NODLo ) *(''! (1*2* 2^M03E^ ) ) * ( M I 1 * 1 * 3^NOUEb ) *|V| (1*5* 54 jNUDE^ ) 
i -M ( 1 * 5»3iNODE5 ) *M ( 1 * 1 * SiNODES ) ) 

Ml ( 1 *4*3ibN0DE3) =MI ( 1 * 3 * 4ii>N0DEb ) - ( M ( 1*2* libNODEb)*M( 1 * 5 ♦ 2^i>NUDEb ) 
b -M ( 1 « 5* 1 ibNODEb ) ii-M (1*2* 2^NODEb ) ) -5^ ( M ( 1 * 1 * 3iN0DEb ) *M (1*4* S^NUDEb ) 
i --IVK 1 *4*3S»N0D£S ) *M( 1 * 1 *5ibNODEb) ) 

M I ( 1 * 4 * 3 SNUDES ) = M 1 ( 1*3* 4i»N0DES ) + ( M ( 1 * 4 * 1 iNODEb ) *M ( 1 * 5 * 2 ^N 0 DEb ) 
i -M ( 1 » 5 * 1 ^NODbb ) «-M ( 1 * 4 * 2^NODEb ) ) * ( Ml 1 » 1 * 3it.|\iUD£b) -«-M ( 1 *2 * Si-NuDEb ) 

s -M( 1 *2 *3‘iiiNOD£S ) *M( 1 * 1 *5iNOD£S) ) 

M I ( 1 * 5*3iNODEb ) = IM ( 1 * 1 * 1 ifcNODES ) *M ( 1 * 2 * 2 itNODEb ) -M( 1*2*1 5>NODES ) * 

S M ( 1 * 1 1 2 ^NOdEo ) ) * ( M ( 1 * 4 * 3ibN0DE-J ) <1*5* 4 bMODEb ) -,v| ( 1 * 5 * 3 SNODE 3 ) 

$ M( 1 *4*4SNOdE3 ) ) 

M I ( 1 * 5 * 3iNODEb ) =M 1(1*3* SbNUUEb ) - I |vi ( 1*1*1 bNUDEb ) (1*4* 2i»NODEb ) 

i -M < 1 * 4 * liNODEb ) *M < 1 * 1 * 2 ibNQOEb ) ) * ( M ( 1 * 2 * 3 ^lN(JUEb ) ■^'^‘■M( 1*5* 4 ^MODEb ) 
£ -M( 1 *5*3^NoD£S) *M( 1 *2*4^N0DES) ) 

Mill *5 *3 *NOd£3 ) = 1 ''! I ( 1*3* S^NOOEb ) + ( M ( 1*1 * 1 ibNOUEb ) -»-M ( 1*5 ♦ 2 SNUDEb ) 
i -M ( 1 * 5 * 1 a^NODES ) *1^1 ( 1 * 1 * 2^N0DES ) ) * ( M (1*2* 3blMODES ) *M ( 1*4 * 4i»N0L)Eb ) 
£ -M< 1 *4*3iii>N0DE5 ) *M( 1 , 2*4£NODEb) ) . 

M I ( 1 * 5 * 3bNOoEb ) =i’'l I ( 1* 3* 5 a>iM JDEb ) + ( M ( 1 * 2 * 1 ^NODEb ) ( 1 * 4 * 2 ^i-NODEb ) 

^ -M < 1 *4 * 1 ^iMODEb ) <1*2* 2^MUD£b ) ) m ( 1 * 1 * 3ibNUDEb ) *M (1*5* 4^NC0£b ) 

S. -M ( 1 *5*3SN0D£S ) ( 1 * 1 * 4iN0DES) ) 

M I ( 1 * 5* 3iNODES ) = M I ( 1 * 3 * SiNODEb ) - ( ( 1 * 2 * 1 iNOQEb ) *i'/l ( 1 * 5 * 2 SNQD£S ) 

S -M < 1 * 5 * l^NODEb ) *M ( 1 * 2 * 2^N0DEb ) ) * < M ( 1 * 1 * 3^N(JDEb ) *M ( 1*4* 4i>NUDES ) 

S -M ( 1 * 4 * 3 ibN 0 DES ) *M ( 1 * 1 * 4 ^i'NODEb ) ) 

M I ( 1 * 5 ♦ 3 itN 0 D£b ) =,''1 1 ( 1 * 3 * S^NODES ) + ( M( 1*4*1 itNODEb ) *ivi ( i * 5 * 2 iN 0 DES ) 
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Table D-3 (Continued) 

t 1 * 5 * 1 i 1 * 4 * 2 ^N 0 DE-::> ) ) * ( M ( 1 * 1 » 3itoNODti3 ) C 1 *2 « 4*NOOE-:j ) 

S -M ( 1 * 2 » 3SNODES ) ( 1 f 1 » 4a>NODES ) ) 

|V| 1 C 1 » I * 4ibN0DE^ ) = UvU 1 * 1 * 2a’N0DE:= ) C i » 2 * 3^NODEi> ) -iM I 1 ♦ 2 * ) * 

a» M ( 1 » i f 3i&N0D£i> ) ) •»■ I M ( 1 *3 * 4^N0DEi3 ) < 1 *5 * 5^-i^NUDE^ ) -M ( I ♦ 5’ 4SislODE^ ) * 

it M{ 1 f3*5itNODES ) ) 

MI ( 1 » 1 *4^b^ 3 dE 6 ) =MI ( 1 * 4* ISNODEia) - ( M( 1 « 1 * 2 SNODEE ) *.vi ( 1 * 3 * 3 SNODES ) 

S -M( 1 *3»2SiMODES) *M ( 1 « 1 i 3 itNODES) )*(M( 1 * 2 • 4i*'N0DEb ) ( 1 *5*5^NoDES) 

* “M ( 1 f5*4itNODi:iS ) *•'"1 ( 1 » 2 » SitNODES ) ) 

|vil ( 1*1 *4itN0DES ) =MI < 1 » 4» 1 ANODES ) + ( M( 1 , 1 * 2itNODEE ) *M ( i * 5 * 3 iN 0 QES ) 

S -M I 1 ♦5*2^N0DE^ ) t 1 ♦ 1 » 3itN0DE^=» ) ) ■«■ ( M ( 1 » 2 » 4itNODEi3) *M ( 1 « 3 » 5--tNODEE ) 
S -M( 1 *3*4it(NOD£S ) *M( 1 * 2 * bitNODEE) ) 

M I ( 1*1 *2^N0DE-^ ) =M I ( 1*4* ia>NODEb ) + C M ( 1 * 2 * EiNuUEE ) ( 1 » 3 * 3itN0D£S ) 

it -iM (1*3* 2itiNODEb ) < 1 » 2 » 3i»NOOES ) ) -Jt- I M ( 1 * 1 * 4i*»iMUL>ES ) *ivi ( 1 * 5 « ) 

* -M ( 1 * 5»4itlMODES ) ( 1 * 1 * 5SNODES ) ) 

M I ( 1 * 1 * 4itNODEo ) =i»l I ( 1 ♦ 4 * 1 ^NUUEE ) - ( M I 1 * 2 » 2 i»NuDEi> ) (1*5* 3itNUDE3 > 

it -M( 1 *5*2itlM0DES ) *M( 1 *2*3itNODES) ) * ( M ( 1 ♦ 1 « 4iwNOUEb ) *M ( 1 *3*5^NuDES) 
a> -M( 1 *3*4itNODt3 ) 1 * 1 * SitoNODES) ) 

M I ( 1 * i * 4itNODE^ ) =M 1 t 1 * 4 * 1 i-NODEb ) + ( M C 1 * 3 * 2^NUDE^ ) ^iv-) <1*5* 3i>NODE3 ) 

S -M ( 1 * 5* 2itNODEb ) ( 1*3* 3iNODE^ ) ) * ( M C 1 * i * 4 :w\JuDE3 ) ( 1 *2 « S^NuQES ) 

S -M ( 1 * 2 ’4S.NODLS ) * 1^1 <1*1* SitNOUEb ) ) 

M I < 1*2* 4itN0D£b ) = ( >'^1 ( 1 * 1 * 1 itNODEb ) *i''l (1*2* 3it>NODE3 ) -M ( 1*2*1 j>NODES ) * 
it M ( 1 * 1 * 3itN0D£sj ) ) * ( M ( 1 *3 * 4itNODE^ ) ( 1 » 5 * ) -H ( 1 ♦ 3 * 4SNUDEi3 ) * 

it )'’j( 1 *3*b^NUDES) ) 

M i ( 1 * 2 ’ 4itNODE'i> ) =/i I ( 1 * 4 • 2itNUDE3 ) - ( i^i ( 1 * 1 * 1 (1*3* 3itN'3DEi=> ) 

it -M ( 1*3*1 ujNQDE-:^ ) (1*1* 3ip)N'JDE^ ) ) ■»■ ( M ( 1*2* 4ipi'^UDEi3 } *)’'i ( 1 * 3 * SipHEEE:^ ) 

-M( 1 *5*4itN0DES ) -«-M( 1 *2*b^N0DE3) ) 

M I ( 1 * 2 * 4itNODE3 ) =1''! I ( 1 « 4 * 2ii>N0uES ) + ( 1 ♦ 1 * 1 i>NU3LE ) ( 1 * 5 » 3i>N0DES ) 

it -M ( 1*5*1 ipiMOOcia } *!'''• (1*1* 3^NuDEii:> ) ) * ( I'^l ( 1 * 2 * 4ii^iMUDE3 ) ^ivi ( 1 * 3 * 5:t>l'JODEi:> ) 

S -M ( 1 * 3*4SNODLS } (1*2* bitNODEb ) ) 

M I { 1 * 2* 4itN0DE3 ) =/j 1 ( 1 * 4 * 2 ipNOOEb ) + ( M ( 1 * 2 * 1 a>lMUOL 3 ) -invi ( 1 * 3 , 3i>iNsjDE iii> ) 

ip -M ( 1*3*1 4>NODL3 ) ( 1 ♦ 2 * 3ibiNUDE3 ) ) -Jt- ( M ( 1 * 1 * 4ii>NoDE^ ) *i«'i (1*5* 5 ^jno je::^ ) 

S -M ( 1 *b*42.i''J0DES ) ( 1 * 1 , SiNODES ) ) 

M I ( 1*2* 4itNODE«3 ) =.^ I ( 1 * 4 * 2it^N0DE3 ) - ( M ( 1*2*1 ipN(JDE5 ) (1*5* 3SNODE5 ) 

it -M ( 1*5*1 itNODEb } ( 1 * 2 * 3ibNODE^ ) ) * ( tM ( 1 ♦ 1 * 4 ^NOuEb ) ( 1 * 3 * 5^:.NUuE^ ) 

S “(VI ( 1 • 3*4a.N0DES ) -»-M ( 1 * 1 * 5^iN<ODES ) ) 

MI ( 1 *2* 4SNODES ) =M I ( 1 * 4 * 2 itNUDEE> ) + ( M ( 1 *3 * 1 itNODEo ) *l'l ( 1 * 5 * 3 a,iMODE 5 ) 

;*> -M( 1 » 5* IS MODE 3 ) *M (1*3* 3*N0DE3 ) ) * ( Ki ( 1 * l * 4it>NuUES ) (1*2* SifiNUUES ) 

S -M ( 1 *2 *4SN0DES ) -«-M (1*1* SitNODEb ) ) 

M I ( 1 * 3 * 4SN0QEb ) = ( M ( 1 * 1 * 1 ibNOUEb ) *M (1*2* 2i»NUDEb ) -M { 1*2*1 itr>JODEb ) * 

it Wi ( 1 * 1 * 2 ^iNUDES ) ) ■«■ ( M ( 1 * 3 * 4 itNUDEb ) *M ( 1 *5 » b^NUDEb ) -Ki ( 1 * 5 * 4 ii>N 3 D£b ) * 

S ivi ( 1 * 3 * 5 ±NOdES ) ) 

M I ( 1 * 3* 44alN(JDE'ii> ) =i"l 1(1*4* 3^f^Oi3Ei::> ) - ( M ( 1 * 1 * 1 it-(NUDEb ) *i-i (1*3 * 2 bNuDEb ) 

S> “M ( 1 * 3 * 1 itNODES ) (1*1* 2 itNODES ) ) * ( |V| ( 1 , 2 * 4bi'^ODEb ) ^I'Al 1 * 5 * S^NUuLb ) 

S “M( 1 *5*4itN0DES) *M( 1 *2»bibN0DES) ) 

MI ( 1 * 3 * 4itN0DE-ii> ) =i''l I ( 1 » 4 * 3 itiM 0 D£b ) + ( M ( 1 * 1 * 1 ioNODEo ) ( 1 * 5 * 2 S 3 N 0 DE -3 ) 
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Table D-3 (Continued) 


5, 1 ,5, liNODtiS) 1 » 1 *2^N0DES) )*tM( 1 » 2*4it=NuDLS)-«-M( 

i “*M ( 1 * 3 « 4*i'JODiiS ) ( 1 * 2 » 5 *I^ 0 QES ) ) 

( 1 ’ 3 ’ 4 j>NODE«J ) I ^ 1 ' 4 ♦ ) + ( |V| ( i * 2 • 1 iuiNUDE'^ ) *!''> ^ 1 


3. 

3> 


S 

a> 

iL 

ifc 

ib 

S 

s 

3. 

ib 

S 

ib 

S 


ib 

ib 

3) 

ib 


“M t 1^3’ l»NODt3 ) I 1 
~M ( 1 • 5»4SN0DES ) *M ( 1 
i''t 1 ( 1 » 3 ’ 4i»NODE’^ ) = I'^l I t 1 
-M ( 1 ♦ 5 » 1 4>N0Dt3 ) { 2 

-M( 1 *3*4ibN0D£S ) *M( 1 

MI ( 1 *3^4snodES ) =;-41 ( 1 

“M ( 1 ♦ 5 ’ 1 ibNODtb ) ( 1 


2 * 2^NOUE^ ) ) * t M ( 1 « 1 » 4^N0UE3 ) -a-M ( 

1 ^ 5i>NODE3 ) ) 

4 * 3^NUDE3 ) - ( M ( 1 , a * 1 ij-NuDb^ ) *r*i ( i 

2 * 2^NuDE3 ) ) * { M ( 1 * 1 * 4^N0DL3 ) ( 

1 , S^NQDES) ) 

4 t 3^(^(^DEi3 ) + (i'^l( 1 *3*1 ibNODt.iiJ ) ■a’M ( 1 
3* 2 ^NUDEi:i ) > * ( M ( 1 ♦ 1 ♦ 4 ^ 1 'MUDEb ) -a-M ( 


-ivi( 1 , 2’4SN0D£S ) ( 1 ♦ 1 ,54>N0DES) ) 

iv'l I ( 1 ♦ 4 » 4ibN0DE^ ) = U’H 1 » 1 » 1 iNODEb ) ^I'^i ( 1 1 2 < 23:>NODEb ) -M ( 1 » 
M ( 1 ♦ 1 » 2 bN 0 DES ) ) * I M ( 1 O » 3bN0DEb ) -K-Ivi ( 1 ♦ E ^ S^NuDEb ) -ivi ( 1 


( 1 f 3 » b 3» N O D E S ) ) 

M I ( 1 * 4 * 4ibN0DE-J J =i4 I ( 1 
“M ( 1 » 3 * 1 ^NODtb ) *M ( 1 
“M ( 1 » 5* 3SN0DES ) -«-M ( 1 
MI ( 1 » 4 »4bNODEb ) =M I ( l 
-M ( 1 * 5 » 1 biMODLb ) *('^1 ( 1 
-M ( 1 < 3 ’3SNO0E5 ) ( 1 

M 1 ( 1 ' 4 ’ 4ibN0DEb ) =i>'l I ( 1 
-M ( 1 » 3 » Ibl-'JODEo ) ^‘•M ( 1 
-M( 1 *5» 33 . nodes ) *M( 1 
MI ( 1 < 4»4SN0dEo ) =M I ( 1 
-M < 1 » 5» 1 ^NODc. 3 ) -a-M ( 1 
-M( 1 OOiNODES) *M( 1 
MIC 1*4’ 4 *i^0d£'^ ) = I*'] I ^ 1 
-M ( 1*5*1 bi\OD ES ) *M ( 1 
-MC 1 »2*3ibN0DES) "a-MK 1 
Ml ( 1 i5*43.NOd£S ) = (MC 1 


4 » 4 ibN(JDEo ) - C Ml ( 1 ♦ I * 1 iwNUDE-iJ ) C 1 ' 
1 » 2iNCiDE^ ) ) * C M C 1 » 2 *3bN0lJEb ) C ; 
2 * SiNODEb ) ) 

4 * 4bNODEb ) + C M C 1 , 1 , 1 ibNODEb ) -a-M C 1 - 

1 » 2^N0DE^ ) ) * C Mi C 1 1 2 * 3*NoDEb ) -a-Mi ( ] 
2* SiNODES ) ) 

4 , 4ifaNODEb ) + C (V| C 1 « 2 * 1 ibNoDL^ ) *M C l < 

2 * 2 bNUDEb ) ) * ( Ml C 1 < 1 » 3it.|MUDEb } *Ml ( ] 

1 *5ibN0DES) ) 

4* 4bNODEb ) - C M < 1 »2* libNODEo ) -a-Mi C 1 , 

2 * 2^m0dE^ ) ) * C |.i C 1 1 1 » 3iijNUljEb ) *Mi C ] 
1 , S^NODEb) ) 

4 * 4 ^NUDEb ) + C Ml C 1 « 3 » 1 d^NOUE^ ) *M1 C 1 < 
3* 2^NODEb ) ) * { M C 1 « 1 » 3bNOUbb ) -a-M ( ] 
1 f SbNODEb ) ) 

1 ^ 1 ibNuDES ) *Mi C 1 » 2 * 2ibNuD£b ) -Mi C 1 » i 


ib Ml C 1 » 1 » 2ibiNjOo E b ) ) ■a- C Ml C 1 » 3 » 3i&N0DEb ; *Mi i 1 ^ o * 4 ifaiNODLii 3 ) — 1*1 C 2 

* M( 1 * 3*4ibi\UDES ) ) I 

M 1 ( 1*5 ♦4i&N0DES ) =i'i i C 1^4* SbNOOEb ) - C Ml C 1 * 1 » 1 ibNODEb ) *M1 C 1 

3> -M CIO* IbNoDES ) -a- Ml ( 10 * 2i-N0DEb ) ) * C M C 1 » 2 * 3it>N0DEb) -a-Ml ( 

a> -MC 1 *5»3ibN0DES ) *M( 1 »2*4ibNODEb) ) 

Ml I { 1 1 5 * 4ibN0DEb ) =i-i 1 C 1 » 4 t 5it.NUDEb ) + C Mi C 1 » 1 * 1 i^NODEo ) -a-Mi C 1 
^ -Ml ( 1*5*1 '"NODcb ) -a- Ml (1*1* 2^''^b'DEO } ) * ( c^i < 1 » 2 * 3ioi^uoEb ) -a-Mi ( 
ib -MCI »3*3ON0 DEo ) ^tMC 1 *2*4ibN0DEb) ) 

M I C 1 * 5 * 4ibNODES ) =ri I C 1 * 4 * 5ibNODEb ) + C Mi C 1*2*1 ibNODEb ) -a-tvi ( 2 

ib -M C 1 * 3 * libNoDES ) *M C 1*2* 2 ibN 0 DES ) ) * C M ( 1 * 1 * JibNOUEb) a-Ml C 

* -MCI * 5* 3bNoQES ) ( 1 * 1 * 4bN0DEb ) ) 

Ml 1 ( 1 * 5 * 4SNODEO ) = ,■'1 I C 1*4* bipNODEb ) - ( )v| C 1 * 2 * 1 ibNODEO > -a-Mi C 1 
ib -M ( 1 * 5* 1 ii^iODEb ) ^Mi (1*2* 2ibNODEb ) ) * C Mi C 1 * 1 * 3ibNuUEb ) * 1^1 C 
s -M c 1 * 3 * 33iNoDEO ) -a-M < 1 * 1 * 4ibNODEib ) ) 


1 * 3* 5*M0DEb ) 

* 3 * abiNiUDEb ) 

L * 5 * Sii’NuoEcb ) 

* 5 * 2 i»NODEb ) 

L * 3 * 53 >NUDEb ) 

•S’E^NODEb ) 

L *2*5bN0DEb) 

E * 1 INODEb ) * 
*5*3S3NUDEb)-«- 

* 3 * 2*M0DE^ ) 

I * 5 * 5 ibl''JuuEb ) 

*5*2ibN0DEb ) 

I * 3 * 5 ^M(jQEb ) 

► 3 * ESNODEi- ) 

[ , 5 * 5 iNODEb ) 

I 5 * 2 £N 0 DEb ) 

L * 3 * SbMUDEb ) 

i5*2ibN0DE^ ) 

L * 2 * SibNODEb ) 

d * 1 SNUDES ) * 

* 5 * 3iNUD£i:5 ) -a- 

» 3 * EiC'NODEb ) 

I *5 *4ibNODEb ) 

* 5*2-bNuDEb ) 

L ♦ 3 * 4^Ni^uEO ) 

. 3*2iN0DEb ) 

, * 5 * 4it>N0DES ) 

I 5 * 2 iNGDEb ) 

. * 3 * 4^NODLO ) 
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Table D-3 (Continiied) 


M I ( 1 « 5 ♦ 4i6N<JDE^ ) =i'''l I < 1 * 4 * S^NODE^ ) + ( M ^ 1 * 3 * 1 ^NODE^ ) ( 1 » 5 * 2 i>NvjDE^ ) 

1 “M ( 1 » 5 ' 1 SsiNOD tS ) ( 1 « 3 ’ 2^f'^EDE^ ) ) ( M ( 1 * 1 * 3^NOQEE ) *i'«'l ( 1 » 2 * 4^NUDE3 ) 

i “M ( 1 » 2 » 3SNODES ) ( 1 » 1 < 4^£>NODES ) ) 

MI ( 1 » 1 * bSNOoEb ) = ( M ( 1 » 1 » 2a-NUDES ) -M-M ( i * 2 » 3ir-rMODEE > -M, I i , 2 »2SN0DEb ) * 

ib M ( 1 1 1 » 3 ^NOdEo ) ) ■«■ 1 1^1 ( 1 < 3 » 4ii>NuDE3 ) -K-M ( 1 ♦ 4 » 5 ii\uUto ) -I'l ( i » 4 » 4iNUD£:3 ) * 

i M ( 1 » 3 » bSiNODES ) ) 

i'’U ( 1 » 1 « S^biNUoEo ) =i^j I ( 1 f 5 » 1 ^NUDE3 ) — ( I''! ( 1 » I * 2 ^f'^*^DEE ) (1*3* 3 ^MuDE^ ) 

i “M I 1 • 3 * 23 jNODE3 ) *M ( 1 * 1 • 3 i£>NUDES ) ) ( Mi ( 1 * 2 * 4^('^UOES ) *i''l ( 1 * 4 « 55 >Nl-DE3 ) 

ib — M ( 1 » 4 * 4S>iMOD ES ) *M ( 1 » 2 « 5i>NUDES ) ) 

I''! I ( 1 ♦ 1 » bSNODE3 ) =.^i I ( 1 » 5 « 1 SNODE3 ) + ( Mi ( 1 ♦ 1 « 2S--NUDE3 ) *|V| ( 1 » 4 , 3 SN 0 DEE ) 

S -M ( 1 « 4 »2a>NODE5 ) ( 1 * 1 » 3i>NUDEb ) ) * ( Ml ( 1*2’ 4it>NUUEb ) * 1^1 (1*3* 5*NUDEb ) 

S “Ml ( 1 *3’4*I^I0DES ) *Ni ( 1*2* SiNODES ) ) 

Ml I ( 1 ♦ 1 » bSNODEb ) =M I ( 1 * b * lit>NOL)ES ) + ( Mi ( 1 * 2 * 2 ^mODE::> ) (1^3, 3il\UDEo ) 

ib — 1 * 3 ’ 2 ^N 0 L)Eb) *Mil 1 * 2 ’ 3it'NUDE3 ) ) * ( Mi ( 1 * 1 * 4 isiMOLL 3 ) C 1 » 4 ♦ 54>NOOE3 ) 

$ -M ( 1 * 4 ’4SNODES ) *M (1*1* S^NODEb ) ) 

M I ( 1 * 1 * ESNOdEo ) =,'U ( 1 ♦ b * 1 bNODES ) - ( Mi ( 1 * 2 * 2iNUDE^ ) <1*4 ♦ 3 S,MUDEb } 

ii> —Ml ( 1 * 4 * 2 ^NoDLO ) (1*2’ 3 ^MuD£b ) ) * ( ( 1 » 1 * 4 i^NUuEb ) -J^ivi (1*3’ b-t'NOuE^ ) 

a. -M ( 1 * 3 ’ 43 ^UD Eb ) (1*1* bi^NUDEb ) ) 

Mil ( 1*1* ba>N JdE-^ ) = ,'1I ( 1 * 5 * liNUDEb ) + ( ,vi ( 1 , 3 * 2'^MoDEb ) *M, (1*4 * 3ii^'»^DEb ) 

ib —Mi (1*4* 2^M0 L)c.o ) ( 1 » 3 » ) ) * ( i'^'j ( 1 * 1 * 4^lNUut:i^ ) *i''i (1*2* ) 

b — M ( 1 « 2 ’ 44 >inI 0DES ) (1*1* bbiNUDEb ) ) 

''’I I ( 1 * 2 * bblNbjjEb ) = ( i'l ( 1 * 1 * 1 ifci'JUbEb ) ■*^•.•1 ( 1 * 2 * 3bi’^o3Eb ) — I'l ( 1 * 2*1 b;\13D£:.C3 ) * 

i M'i(l*l*3bNODE'-’))*( *''' ( 1 *3*4bNuDEb)*i''i( 1 *4*3-^MbDEb)— Mi( 1 »4»4Sinj oDE^)* 
b Ml ( 1 , 3 * bbNUoES ) ) 

Ml I ( 1 * 2 * b4DNODEb ) =ri I ( 1 * b* 2 bNOUEb ) - ( iv| ( 1 * 1*1 ibNODEb ) *mi (1*3 * 3 ibiMUDEb ) 

b -Ml ( 1 * 3 ’ 1 bNQDEo ) *Mi ( 1 * 1 * 3 bNODEb ))-«-( Mi ( 1 * 2 ’ 4bNuUEb ) (1*4* S^MUDEb ) 

i -M { 1 , 4 » 4b nodes ) (1*2* bbNODES ) ) 

Ml I ( 1 * 2 * bbNUDEb ) =.-'i I ( 1 * b * 2 bNODEb ) + ( Mi ( 1 ♦ 1 * l i.NUDEb ) ^iv, ( 1 * 4 » 3 bNuDE :3 ) 

1 -M ( 1*4*1 ajNODEb ) ^I'^i (1*1* 3^NUDEb ) ) * ( i''l ( 1 ♦ 2 » 44JiNjODEb ) *i-i (1*3’ SbNODEb ) 

b -M ( 1 * 3 *4bNQDES ) -«-Ml ( 1*2* SbNODEb ) ) 

Ml I ( 1*2’ bbNUDE-3 ) = i''i I ( 1 ♦ b » 2bNOL)Eb ) + ( I'^i ( 1*2*1 bi\UDEb ) *i''l (1*3* 3bNvJLEC3 ) 
b -Ml ( 1*3*1 ibNODt-b ) -«-Mi (1*2* 3bNOOEb ) ) ( Mi ( 1 * 1 ♦ 4^NODEb ) *M (1*4* S^MlDuc^b ) 

i -,-71 (1*4 * 4bNODES ) *Mi (1*1* SbNODEb ) ) 

Ml I ( 1 * 2 * bbNOoES ) =m 1 I ( 1 * 5 * 2bNODEb ) - ( Ml ( 1 * 2 * 1 bNUDEb ) *M (1*4* 3bi\CDEb ) 

b — M ( 1 * 4 ’ 1 bNODLb ) *Ml ( 1 * 2 ’ 3bNOQEb ) ) * ( Ml ( 1 * 1 * 4bNODEb> ) (1*3* Sbi'-lODEb ) 

b -M( 1 *3*4bN0DES ) *M ( 1 * 1 * SbNODES) ) 

M I ( 1 * 2 * bbNUDEb ) =Mi 1 ( 1 * b * 2bNUDEb ) + ( Mi ( 1 » 3 * 1 iciNODEb ) (1*4* 3 bNUDEb ) 

b — ivi ( 1 , 4 , 1 bNODEb ) ■> 4 'Mi ( 1 * 3 * 3 bN 0 DEb ) ) * ( Mi ( 1 * 1 * 4 bNUDEb ) *Mi ( 1 * 2 * SbNODEb ) 

b — M ( 1 * 2 * 4 bNOD tS ) ■^Mi (1*1* SbNODEb ) ) 

Ml I ( 1 * 3 * SSNOdES ) = ( M ( 1 * 1 * 1 bNUDEb ) *Mi ( 1 * 2 * 2 bNODES ) -Ml ( 1 * 2*1 bNODES ) * 

b .'' 1 ( 1 * 1 * 2bNODEb ) ) * ( N ( 1 * 3 * 4bNUDEb ) *M 1 ( 1 * 4 * SbNUDEb ) -M ( 1 * 4 * 4blM0DEb ) * 

S M { 1 * 3 * bbNODES ) ) 

Ml I ( 1 ♦ 3 * SSNOdES ) ^mI 1 ( 1 * b* 3bNODEb ) - ( M'( 1 * 1 * 1 bNUDEb ) (1*3 * 2bNUDE^ ) 

b — M)( 1 * 3 ’ IbNODuS) ( 1 * 1 * 2 bl'Jv_;DEb ) )*(ivi( 1 * 2 * 4 bNUDLb)-^^-Ml( 1 * 4 * S^telNODEu ) 

b —Ml (1*4 * 4bNODES ) (1*2* SbNODEb ) > 

Mil ( 1 *3* SbNODEb ) = MI ( 1 * 5 * 3bNOOEb) + ( Ml ( 1 * 1 * 1 bNUDEb ) *iv| ( 1 , 4,2bN0DEb ) 

(Continued) 


169 
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s -M < 1 » 4 » ia>NODbS ) ( 1 » 1 * a^bNODES ) ) ■»■ ( M ( 1 » 2 » A^bNOUES ) *fvi ( 

b -MCI O»4bi\0DES) ^M( 1 »2»5bN0DES) ) 


M I ( 1 

it ( 
S “M( 


^3’5b^ OoEi ) =M I (■ I 


5 5 3bNQDES ) + ( M ( 1 « 2 » 1 bNUDEb ) C 1 


M3* lb<MCDE5) 1 »2*2bN0DEcj) ) * ( iv/i ( 1 * 1 ♦ AbNODES ) *|V| ( 
1 , A*ASN0DES ) ^tivi ( 1 , 1 * SiNODES) ) 


( 1 i3*E-pN0dE^ ) =.^il ( 1 
i -,M ( M A * 1 bl'JODEcs } *1''! ( 1 

b -M ( 1 , 3 « AbIsjODLS ) *M ( 1 

I ( 1 ♦ 3 »KbNODEi) ) =i"i i ( 1 

iip -[71 ( M A * 1 bNQL) C 1 

b -M ( M 2 * AbNOOcS ) ( 1 

M I ( 1 » A * SbNUDEo ) = ( M ( 1 

b M ( 1 ♦ 1 * 2 ^NOdE-3 ) ) ^ M ^ 

b ( 1 » 3 * EbNOD E S ) ) 

M I ( M A * EbNOpEb ) =i''i I ( 1 
b -M ( 1 * 3 » 1 bNODES ) ( 1 
b -M. ( 1 , A * 3bN0DES ) ( 1 

MI ( 1 » 4 »bbNODEo ) =MI C i 
i, -M ( 1 ♦ A * 1 bNQO ES ) ■>'^1''! ( 1 
b -M( M3*3bi-xiODES ) ’'‘■i-K 1 
i''i I ( 1 ♦ A * b^i'^Oot A ) = i-l 1 ( 1 

b -|V| ( 10*1 ^NOD lS ) «-M ( 1 

b ~1V| (1,4 ♦ abiMODES ) "'^1''' ( 1 

M I ( 1 * A * bbNUDEb ) =M I ( 1 
S -M { 1 * A * MNODlS ) C 1 
b -M( 1 * 3*3bi\i0DES ) ^i'l ( 1 
MI ( 1 * A * bSN OD E ) = Ml ( 1 
S -M ( 1 ♦ A » 1 bl'JODES ) 1 

b -M { 1 * 2 * 3bNODES ) *M ( 1 
Ml { 1 ♦ 5 * bbNuDEo ) = ( M ( 1 
b I''! { M 1 * 2-bNbDE-> ) ) "^ 1 f’’! ( 
b K1 { 1 O * AblNiODES ) ) 

M 1(1*5* bbNODEE ) =)''i I ( 1 
b “M ( 1 * 3 * 1 biMODES ) ( 1 

b -M( 1 ,A*3bN0DES) 1 
M I ( 1 * 5* bbNUDE^ ) =M I ( 1 

b -M ( 1 * A * 1 bNOptb ) C 1 

b ~M ( 1 * 3 * 3bN0DES ) ■’'''■M ( 1 

M I ( 1 * 5* SbNODEE ) I ( 1 

^ -M ( 1 * 3 ♦ 1 bNODEE ) ■’"‘•M ( 1 

b -M ( 1 * A *3bNODES ) ( 1 

1^1 I ( 1 *5*bb^ 3 dEo ) =i‘'ll ( 1 

b ~M ( 1 * A » 1 ai'MOD Lo ) ( 1 

b -M( 1 *3*3biN0DES ) *M ( 1 

1'^) I ( 1 * 5 * b^N OdE C:> ) = M I ( 1 


* 5* 3bNODEb ) - ( M ( 1 ♦ 2 * 1 bNODEb ) ( 1 

* 2 * 2bl'^ODEb ))'W’(i7|(Mi» AbNObES ) *M ( 

* 1 * SbNODES ) ) 

* 5 * 3bNUDES ) + ( I'l ( 1*3*1 bNUDEb ) ( i 

* 3 * 2biNiOiJE'i^ ) ) * ( i''i ( 1 * 1 * Abl'-g CJbEb ) ( 

* 1 * SbNODEb ) ) 

*1*1 bNODEb ) -X-|v| (1*2* 2bNUDEb ) “>^1 ( 1 » 
1*3* 3bN0DEb ) ^-M ( 1 ♦ a * SbNODES ) C i 


* b * AbNGUEb ) ~ ( I'-'t ( 1*1*1 bNODEb ) ’^1'“': ( 1 

* 1 ♦ 2bi\lODE3 ) ) * { M ( 1 » 2 ’ 3bNODEb ) -Ji-M ( 

* 2* bbNODES ) ) 

* 5 , AbNODES ) + ( M ( 1 * 1 * i ^u^NODEb ) - i^i ( 1 

, 1 , 2bfMoDEb ) ) * ( M ( 1 * 2 * 3bi\UDEb ) ( 

*2* SbNODES) ) 

* ^ * AbNODE b ) + ( I'^i ( 1*2*1 bNUDEb ) ( i 

* 2 * 2bNODEb ))-«■( 1^1 ( 1 * 1 * 3bNODES ) ( 

* 1 * SbNODES ) ) 

* D * AbNODES ) - ( M (1*2*1 bNoDEb ) *i/| ( i 

* 2 ’ 2bNODEb ) ) * ( i^i ( 1 * 1 * SblNOUES ) ( 

, 1 , SbNODES ) ) 

* b* AbNODES ) + ( M ( 1 * 3 * 1 bNODES ) -x-M ( 1 

* 3 * 2bN0D£S ) ) * ( M ( 1 * 1 * 3 bNODES ) '--2. ( 

* 1 , SbNODES ) ) 

*1*1 bNODES ) "jY iv| (1*2* 2bN0UES ) — l-i ( 1 * 
1*3* SbNODES ) -«-|vi (1*4* 4:j:;iNODEb ) -|V| ( i 

* b * bbNuDES ) - ( M ( 1 * 1 * 1 ANODES ) ( i 

* 1 * 2biMODES ) ) * ( M ( 1 * 2 * SbNODES ) *M ( 

, 2* AbNODES ) ) 

* 5 * SbNODES ) + ( 1^1 ( 1 * 1 * 1 bNODES ) -^|vt ( i 

* 1 * 2bNODES ) ) * ( M ( 1 * 2 * SbNODES ) *M ( 

, 2* AbNODES ) ) 

* b* SbNODES ) + ( M ( 1*2*1 bNODES ) *M ( 1 

* 2 * 2bN«-^DES ) ) * ( M ( 1 * 1 ♦ SbNODES ) ( 

, 1 * AbNODES ) ) 

« 5 * SbNODES ) ~ ( I''! { 1*2*1 blMODES ) *ivi ( i 

* 2 * 2bNODES ) ) ^ ( M ( 1 * 1 * SbNODES ) -K-,v, ( 

* 1 * AbNODES ) ) 

* 5* SbNODES ) + ( N ( 1 * 3 * I bNODES ) *M ( 1 


1 * 3 * SbNODES ) 

* 3*2bN0DES) 
1*4* SbNODES ) 

* A * 2bNODE :i3 ) 

1 * 3 * S^NUDES ) 

* A *2b NODES ) 
1*2* SbNODES ) 

2 * 1 bNODES ) 

*A*3bNOD£-;^)* 

* 3 * 2bN0DES ) 

1 *A*5bN0DEs) 

,4*2bNOUES) 

1 *3*5^ NODES) 

* 3 * 2bNUUES ) 

1 *A* SbNODES) 

* A * SbNODES ) 

1 * 3 * SbNODES ) 

* A * sbNOuES ) 

1*2* SbNODES ) 

2 » 1 bNODES ) 

* A * sbNoDEs ) * 

* 3*2bN0DEs ) 

1 * 4 * A^NODEb ) 

, A * SbNODES ) 

1 *3*A''^NUDEi^) 

* 3 * SbNODES ) 

1 * A * AbNODES ) 

* A , SbNODEb ) 

1 * 3 * AbNODES ) 

* A * sbiNiODEs ) 
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Table D-3 (Continued) 

-M ( 1 • 4 » laiNODtS ). t 1 * 3» 2ii»NGUt.S ) ) * ( M ( 1 * 1 » 3*i'*ODE3 ) ^-ivi ( 1 » A^NOGEb ) 

* -M ( 1 » 2 * 3^N0D£S ) ( I ♦ 1 » A^NODES ) ) 

DET ( 1 SNUUEb ) =i'H 1 * 1 « 1 ^i''JODES ) *M I ( 1 « 1 » 1 SNODES ) 

■f ri ( 1 » 1 » 23iNOUES ) *M I ( 1 ♦ 3 » 1 SNODES ) 

+ N ( 1 > 1 * 3SNODES ) -«-M I ( 1 « 3 < 1 SFnIODES ) 

4 1 +,v| ( 1 , 1 , 4ibiMODES )*M i ( 1 » 1 » ASNODES ) 

+ 1''! { 1 » 1 * 5SNODES ) *M I ( 1 ♦ 5 * 1 SNODES ) 

NMi = nodes- 1 
NM2=N0DES-2 
DO 50 I = 1 « 5 
DO 50 J=l»5 


i‘1 1 ( 1 « 1 ♦ oysi'jODEo ) = 1-1 1 ( 1 « I » Denude b > /det ( i si\odes ) 

50 continue 
I TEN=0 

10 continue 

I TeN= I TEN+ 1 

DO 100 I = 1 » 5 

uP( 1 » I ) = ( 1 ) *u ( 1 * I ) 

UP ( 1 « I ) =UP ( 1 9 I ) I ( 1 9 I 9 1 ) * ( In ( 2 ’ 1 * 1 ) *u ( 2 ’ 1 ) +IN ( 2 » 1 <• 2 ) *U ( 2* 2 ) 

S +N < 2 ’ 1 » 3 ) *U ( 2 » 3 ) +N <2M»A)*U(2*A)+N{2«l95)*U(2»5) ) 

UP ( 1 9 I ) =UP ( 1 > I ) - I( l9 I*2)*(iNi(2’2*l) ■^^U ( 2 » 1 > +N (29292)*U(2»2) 

S AN ( 2 9 2 ♦ 3 ) ^U ( 2 ♦ 3 ) +N ( 2 ’ 2 ♦ 4 ) *U ( 2 9 A ) -fN ( 2 * 2 9 5 ) -^^-U { 2 » 5 ) ) 

UP (l9l)=UP(i9l) I ( 1 9 I 9 3 ) * ( N ( 2 9 3 9 1 ) -«-u (2*1) AN (29392) *U(2’2) 

S AN ( 2 9 3 9 3 ) *U ( 2 9 3 ) AN ( 2 9 3 9 A ) -«-U ( 2 9 A ) AN (29395) *0(295) ) 

UP ( 1 9 I ) =UP ( 1 9 I ) - ,J*|V| I ( 1 9 I 9 4 ) * ( N ( 2 9 A 9 1 ) *U ( 2 9 1 ) AN (29492) *U(292) 

S AN ( 2 * A 9 3 ) *U ( 2 * 3 ) AN ( 2 94 9 A ) *u ( 2 9 4 ) an ( 2 94 9 5 ) *0 ( 2 9 5 ) ) 

UP ( 1 9 I ) =UP ( 1 9 I ) - I ( i 9 I 9 5 ) * ( N ( 2 9 5 9 1 ) *u ( 2 9 1 ) AN ( 29 59 2 ) *u ( 2’ 2 ) 

S AN ( 2 9 59 3 ) ^U ( 2 9 3 ) an ( 2 9 5 9 4 ) AU ( 2 9 4 ) AN ( 2 95 9 5 ) *U ( 2 9 5 ) ) 

UP ( 1 9 I ) =UP ( 1 9 i ) A iij A ( |V‘ I ( 1 9 I 9 1 ) AD ( 1 9 1 ) AM I ( 1 9 I 9 2 ) *D ( 1 9 2 ) A 

S I ( 1 9 I 93 ) *D ( 1 , 3 ) A.M I ( 1 9 I 9 4 ) *D ( 1 9 4 ) AM I ( 1 9 I 9 5 ) *D( 1 9 5 ) ) 

UP ( 2 9 I SNM2 ) = ( 1 .-vV ) *U ( 29 I SNM2 ) 

UP (291 SNi''l2 ) =uP ( 2 9 I ^NM2 ) — »/wAivi I ( 2 9 I 9 1 4jNM2 ) * ( L ( 1 9 1 9 1 UiNi'/i 2 ) aij ( 1 9 1 Givi'>';2 ) 
s AL ( i * 1 9 2SN|Vi2 ) Au ( 1 9 2SnM2 ) AL ( 1 ♦ 1 9 3 SNM 2 ) *u { 1 9 3 SNM 2 ) AL ( 1 9 1 9 4ipN|vi2 ) A 
S U ( 1 9 4SNi^i2 ) aL ( 1 9 i 9 5SNN2 ) *U ( 1 9 5ANi^i2 ) ) 

UP ( 2 9 I SNM2 ) =UP (291 SNM2 ) - WAN I ( 2 9 I 9 1 *NM2 ) A ( n ( 3 9 1 9 1 u^NN2 ) *0 ( 3 9 1 SNM2 ) 
s AN (3919 2^Nm2 ) Au ( 3 9 2SNM2 ) AN ( 3 * 1 9 3SNM2 ) *U ( 3 9 3 SNM 2 ) AN ( 3 9 1 9 4-pNiVi2 ) * 
S 0(39 4SNM2 ) aN ( 3 > 1 9 5SNM2 ) AU ( 3 9 5SNiH2 ) ) 

UP ( 2 9 I SNM2 ) =UF ( 2 > I a>NM2 ) - WAp, I ( 2 * I 9 2i»NM2 ) * ( L ( 1 9 2 9 1 -i.NM2 ) *0 ( 1 9 1 SNM2 ) 

S AL < 1 9 2 9 2SNM2 ) AU { 1 9 2SNM2 ) AL ( 1 9 2 9 3^NM2 ) *u ( 1 9 3u.N|V|2 ) AL ( 1 9 2 9 4SNM2 ) A 

S U ( 1 9 ASNM2 ) aL ( 1 » 2 9 5SNi^)2 ) *U ( 1 9 5SNM2 ) ) 

UP ( 2 9 I SNM2 ) =UP ( 2 » I SNM2 ) - wAM I ( 2 9 I 9 2 SN|vi 2 ) a ( n ( 3 9 2 9 1 ^jNi''I2 ) *0 ( 3 9 1 SNM 2 ) 
S AN ( 3 9 2 9 2ibNM2 ) AU ( 3 9 2SNM2 ) AN (3929 3SNM2 ) *U ( 3 9 3a-NM2 ) AN (3939 4a>NM2 ) * 
S 0(39 4SNM2 ) aN ( 3 ♦ 2 9 5SNM2 ) AU ( 3 9 5SNM2 ) ) 

UP ( 2 9 I SNPI 2 ) =uP ( 2 » I SNM2 > - wan I ( 2 * 1 9 343 NM 2 ) A- (L ( 1 9391 SNP 12 ) AU (Ml it.NM2 ) 
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Table D-3 (Continued) 


+L t 1 » 3 » 2i»NiVi2 ) *u t 1 ’ 2ibiN|Vi2 ) + 1 - t 1 » 3 * 3 i>NM 2 ) *0 ( l ♦ 3*i\iiVi2 ) +L ( 1 * 3 * 44 >i\M 2 ) * 

£ U ( 1 » 4i^N|vi2 ) +L ( 1 ♦ 3 * 5it>NM2 ) *U ( 1 » 5i»N^i2 ) ) 

'^P (2*1 5>Ni''12 ) =UP ( 2 » I ^Niv|2 ) “ W*|V| I < 2 * I * 3^Nlvl2 ) * ( N ( 3 * 3 » 1 ij3NM2 ) *U ( 3 * 1 $NM2 ) 
i> +N ( 3 • 3 *2^NM2 ) < 3 * 2ilMM2 ) +N( 3*3* 3^Nivi2 ) *U ( 3 * 3^NM2 ) +1^ (3*3 « 4it>NKi2 ) * 

* U ( 3 ♦ 4ifcNM2 ) +jNi ( 3 > 3 * 5*NM2 ) -«-U ( 3 * b^NM2 ) ) 

UP (2*1 3jiNti''l2 ) =<jP ( 2 » I it>NiVl2 ) 1(2*1* 4^i'^M2 )'^(t_( 1 «4» 1 iNI''i2 ) *U ( 1 ♦ 1 ii»Nivi2 ) 

+L< 1 *4*2^6N!V12) *u( 1 * 2il.N^•,2 ) +L ( 1 *4 » 3^Nlv'i2 ) -»-U ( 1 ♦ 3^»INM2 ) +U ( 1 .4*4 *Njv12 )* 

4 U { 1 » 44NM2 ) + L ( 1 ♦ 4 * 54'Nr''l2 ) *U ( 1 « 54NH2 ) ) 

UP (2*1 S>i'^i'^t2 ) =up ( 2 » I ^.MM2 ) -W*,V| I ( 2 * I * 4i>Nl''i2 ) * ( N ( 3 » 4 » 1 4jN^ 2 ) *u ( 3*1 4NM2 ) 

* +N ( 3 * 4 * 2a^NM2 ) •"'‘•U ( 3 * 2iNM2 ) +N ( 3 * 4 * 3^NM2 ) *U ( 3 » 34NM2 ) +N ( 3 « 4 » 44 NM 2 ) ^ 

^ U ( 3 « 4Si\yi2 ) +N ( 3 > 4 * 5^NM2 ) ->^U ( 3 » 54NiVi2 ) ) 

up (2*1 4 Ni'^2 ) = '-'P ( 2 » I 4iN|Vi 2 ) — 1(2*1* 5^Nl''l2 ) * ( L ( 1 * 5 * 1 ‘-biNi|V|2 ) ^*U ( 1 ♦ 1 4N|vi2 ) 

^ +L ( 1 * 5 * 2^i'^M2 ) ( 1 * 2^''^i'’i2 ) +1- (1*5* 3^i^i''i2 ) "^U ( 1 ♦ 3 :t.i\ivi 2 ) -hL_ (1*5 ♦ 4•UJ^^Jl'l2 ) * 

4 U ( 1 » 44NM2 ) +L ( 1 ♦ 3 * 5iNiVi2 ) *U ( 1 » 5iNM2 ) ) 

^P (2*1 iiNi'i^i2 ) =UP ( 2 > I ^Niv' 12 ) - 1(2*1* 5^iNi'^i2 ) "^ ( N ( 3 * 5 * 1 -t-Niv'i2 ) *U ( 3 * 1 ^fcNivi2 ) 

i +N ( 3 * 5 * 2iNM2 ) ( 3 » 2^ImM 2 ) +N (3*5* 34fMi''i2 ) *U ( 3 * 3^iNwl2 ) +N ( 3 » 5 ♦ 4^Nwi2 ) * 

4 U ( 3 f 4ilMf’^2 ) +fM ( 3 * 5 « 5iNjvi2 ) *U ( 3 * 54NM2 ) ) 

^P (2*1 4Niv( 2 ) =UP (2*1 ^Nr'>'i 2 ) +i\>* ( iw I ( 2 * I * 1 4 Niv| 2 ) (2*1 4Ni''l2 ) +M 1(2*1* 2 ^Nivi 2 ) 

4 ( 2 * 2-°i''J^i2 ) +i'l 1(2*1 * 34NM2 ) *D ( 2 * 3i'NM2 ) +Ki 1(2*1’ 44NIV-I2 ) *D ( 2 • 4 a>NM 2 ) 

4 +,vi I ( 2 * I * 54NM2 ) *D ( 2 *54NM2 ) ) 

UP( NODE^* I ) = ( 1 .-vJ) *U( NODES* I ) 

UP( NODES* 1 ) =UP { NODES* I (NODEi^* 1*1) * ( L ( NM 1 * 1 * 1 ) ^'^U ( NM 1 * 1 ) 

4 +L ( NM 1 * 1 , 2 ) * O ( NPi 1 * 2 ) +L ( NM 1 * 1 * 3 ) *U ( Ni^l 1 , 3 ) +L ( N|V| 1 *1*4) 

4 -itj ( nM 1 * 4 ) +L I I 1 * 1 * 5 ) *U ( NM 1*5)) 

uP ( NUDL^ * I ) =UP ( I'^ODEb * I ) — I ( NoDEb * 1 ♦ 2 ) "^ ( L ( Ni<'i 1 * 2 * 1 ) *U ( Ni'^i 1 * 1 ) 

4 +L(NM1 * 2 * 2 ) ^'^ ( NM 1 * 2 ) +L ( NM 1 ♦ 2 * 3 ) *U ( NM 1 * 3 ) +E ( NM 1 *2*4) 

4 *u ( NM 1 * 4 ) +L ( i^M 1 * 2 * 5 ) *U ( NM 1*5)) 

oP ( NODE-^ * I ) =oP ( N JDE^ * I ) “* w*ivi I ( NUDE^ * I * 3 ) * ( L ( 1 * 3 * 1 ) '^U ( NiM 1*1) 

4 +L ( NM 1 * 3 * 2 ) '^ U ( Ni’i 1 * 2 ) +E ( NiVi 1 * 3 * 3 ) *U ( Ni''i 1*3) +L ( iNi-i 1 *3*4) 

4 ( NM 1*4) +L ( ■ 'li''! 1 *3*5) *U ( NM 1*5)) 

UP ( NODES * I ) =UP { IM ODES * I ) -w-«-M I ( i\UDE^ * I *4 ) * ( L ( NM 1 * 4 * 1 ) *U ( NM 1 * 1 ) 

4 +L ( NM 1 * 4 * 2 ) ( Ni^'i 1 * 2 ) +E { NMi 1 * 4 * 3 ) *U ( Ni^i 1 * 3 ) +L ( NM 1 * 4 * 4 ) 

4 *U( NMl *4 ) +L ( imMI * 4 * 5 )'-5<-U ( NM 1 *5 ) ) 

UP ( NODES * I ) = UP ( /iODES * 1 ) -W*M I ( NODES * 1 * 5 ) * ( L ( N^i 1 *5*1) *U ( NM 1*1) 

4 +L ( N‘^l 1 * 5 ♦ 2 ) ( Ni'’i 1*2) +E ( NM 1 *5*3) *U ( NM 1*3) +L ( nMI 1 *5*4) 

4 *U ( Nl*l 1*4) +L ( 1 4‘'1 1 *5*5) *U ( NMi 1*5)) 

UP ( N0D£-:3 * I ) =UP ( I'^UDLS * I ) + kj* ( Ml I ( NODES *1*1) *D ( NuDES * 1 ) +M I ( NODES *1*2) 

4 *D ( NOUc_S « 2 ) +i'i I ( imUDES *1*3) *D ( NUDE^^a * 3 ) +M I ( NODE^ *1*4) '^D ( NODE-^j * 4 ) 

4 +M 1 ( nudes * I > 5 ) *D ( NODES * 5 ) ) 
luo continue 

DU 350 1=1*5 

DIF( 1 * I*NODeS) = JP( 1 * I4N0DES)-U( 1 * I 4NODES ) 
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Table D-3 (Concluded) 


3b0 continue 

RMS=0 • 

DEL=0. 

DO 360 1=1*5 
DO 360 J=1 *nODE^ 

DEL = DEL+D1 F ( J * I ) *0 1 F ( J * I ) 

RmS = RMS+UP ( J > I ) -«-UP ( J * I ) 

360 continue 

DEL=S0RT ( DEL ) 

RMS = SORT ( RMS ) 

test=del/rms 

WRITE(6»4£ D)UP, J»DIF 
450 FORMAT! 9E1 0. 3 ) 

W/Rl TE ( 6* 500 ) I TtR * test » DEL 1 RMS 
500 format ( I 5* 3e10 • 3 ) 

DO 400 1=1*5 

U ( 1 * I iNODES ) =UP ( 1 « I ANODES ) 

400 CONTINUE 

I F ( TEST - LE . ePv:> • OR . I TER. GE . MAX 1 ) RclTUKN 
GO TO lU 

end 
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